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Abstract 


NORTH,  JAN  A.,  Spatial  Information  Content  Analysis  of  Optical  Polarimetric 
Imagery.  Typed  and  bound  thesis,  176  pages,  11  tables,  113  figures,  and  123  equations. 
May  1995. 

In  this  study  the  spatial  information  content  of  optical  polarimetric  (Stokes  vector)  imagery 
of  both  synthetic  and  natural  scenes  is  analyzed.  The  intent  of  this  study  is  to  evaluate  the 
utility  of  polarimetric  imagery  in  the  extraction  of  spatial  information  content  from  scenes 
containing  spatially  correlated  features  that  can  be  described  by  their  spatial  (two- 
dimensional)  power  spectra. 

This  study  was  conducted  in  four  phases;  1)  construct  an  imaging  polarimeter  system  for 
the  purpose  of  collecting  representative  imagery  containing  both  synthetic  and  natural 
scenes  of  spatially  correlated  features  illuminated  under  clear  sky  conditions;  2)  process 
the  collected  image  data  to  create  the  Stokes  vector  images  and  derivative  images;  3) 
construct  a  polarimetric  imaging  simulation  for  the  purpose  of  comparing  synthetic  data 
with  actual  results;  and  4)  analyze  both  the  synthetic  and  actual  Stokes  imagery  for  the 
purpose  of  evaluating  the  spatial  information  content  of  the  processed  scenes. 

Analysis  of  the  study  results  produced  five  main  conclusions:  1)  characteristic  curve  co¬ 
calibration  of  the  input  images  is  important  for  the  accurate  calculation  of  the  Stokes 
parameters  and  derivatives;  2)  the  polarimetric  difference  (D)  image  provides  a  direct 
measurement  of  sensitometric  co-calibration  error;  3)  the  polarization  ellipse  orientation 
angle  (T)  image  histogram  can  provide  an  indirect  measure  of  the  azimuthal  difference 
between  the  principal  plane  and  the  polarizer  reference  plane;  4)  the  percent  polarization 
(P)  and  orientation  angle  (T)  components  provide  more  stable  estimates  of  spatial  power 
spectral  density  (PSD)  compared  with  the  unpolarized  intensity  (I)  component;  and  5)  the 
PSD  stablities  of  the  P  and  T  components  demonstrate  a  relative  insensitivity  to  imaging 
geometry  compared  with  the  I  component. 

The  results  of  this  study  provide  a  refined  understanding  of  the  potential  contributions  that 
polarimetric  imagery  can  provide  in  the  analysis  of  spatial  information  content  from  optical 
imagery. 
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1.0  Introduction 

1.1  Purpose 

The  purpose  of  this  study  is  to  obtain  an  improved  understanding  of 
the  contribution  that  sensed  polarization  can  provide  in  the  analysis  of 
spatial  information  content  from  optical  imagery.  Toward  that  end,  the 
primary  focus  is  to  obtain  both  qualitative  and  quantitative  evaluations  of 
the  utility  of  optical  polarimetric  (Stokes  vector)  imagery  in  the  extraction 
of  spatial  information  content  from  scenes  containing  spatially  correlated 
features  that  can  be  described  by  their  spatial  power  spectra.  Polarimetric 
imagery  containing  both  natural  and  synthetic  scene  features  are 
analyzed. 

1.2  Objectives 

This  objectives  of  this  investigation  are  described  in  four  phases: 

Phase  1  -  Image  Data  Collection 

Construct  an  imaging  polarimeter  system  using  a  4-lens  35-mm 
camera  with  linearly  polarizing  filters  differentially  rotated  about  their  optic 
axes.  Execute  a  controlled  experiment;  collect  representative  polarized 
imagery  containing  both  synthetic  and  natural  scenes  of  spatially 
correlated  features  illuminated  under  naturally  polarized  (clear  sky) 
illumination  conditions. 

Phase  2  -  Image  Data  Processing 
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Process  the  differentially  polarized  image  sets  to  create  calibrated 
Stokes  parameter  images  and  Stokes  derivative  images;  convert  the 
photographic  images  into  digital  images;  co-register  the  digital  images;  co¬ 
calibrate  the  images  for  relative  sensitometry;  and  then  mathematically 
combine  the  co-registered,  co-calibrated  digital  images  to  create  the 
Stokes  parameter  and  Stokes  derivative  imagery. 

Phase  3  -  Image  Data  Simulation 

Synthesize  a  reasonable  imaging  model  which  can  simulate  the 
effects  of  resolving  differentially  polarized  imagery  of  a  mirror-reflected 
skydome  under  various  imaging  geometries.  Exercise  the  model  to  create 
synthetic  images  which  can  be  compared  with  actual  calibration  images  of 
the  mirror-reflected  skydome. 

Phase  4  -  Image  Data  Analysis 

Analyze  the  errors  associated  with  calibrating  the  Stokes 
parameter  and  Stokes  derivative  imagery  from  digitally  processed,  4-lens 
35-mm  polarized  photographs. 

Analyze  and  compare  both  the  synthetic  and  actual  polarized  image 
features;  specify  sub-images  within  each  image  that  contain  similar  in¬ 
scene  content  and  texture;  transform  these  sub-images  into  their  spatial 
power  spectra;  analyze  and  compare  the  resulting  spectra  to  determine 
the  effects  that  imaging  geometry  may  have  on  the  stability  of  spatial 
spectral  estimation  for  each  of  the  polarized  image  components  (i.e.,  the 
Stokes  parameter  and  Stokes  derivative  images). 
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1 .3  Motivation 


The  prime  motivation  for  this  study  is  to  utilize  a  large  body  of 
current  research  of  one  spatially  correlated  natural  feature,  v\rind-driven 
waterwave  surfaces,  as  the  basis  for  the  study  of  other  natural 
background  and  clutter  features,  definable  by  their  spatial  power  (Wiener) 
spectra,  under  polarized  optical  imaging  conditions.  A  secondary 
motivation  is  to  consider  the  potential  utility  of  polarized  waterwave 
scenes  as  natural  opportunistic  targets  for  synoptic  image  calibration.  This 
motivation  is  based  on  several  unique  aspects  of  water  and  waterwaves 
enumerated  below.  A  final  motivation  is  to  demonstrate  the  feasibility  of 
extending  this  analysis  to  more  sophisticated  imaging  systems,  using  the 
relatively  simple  techniques  employed  in  this  study. 

There  are  several  reasons  for  selecting  natural  waterwave  surfaces 
(and  their  simulcra)  as  the  initial  object  of  polarized  image  analysis: 

1)  Terrestrial  surface  water  is  ubiquitous;  fully  three  quarters  of  the 
Earth's  surface  is  covered  with  water.  For  airborne  and  spaceborne 
sensors,  this  represents  a  large  window  of  opportunity  to  image  this 
natural  feature,  even  in  the  presence  of  clouds. 

2)  In  comparison  with  other  terrestrial  surface  features,  water 
surfaces  are  optically  homogeneous:  only  two  radiance  functions,  one  for 
surface  reflection  and  one  for  upwelling  subsurface  refraction  (i.e.,  volume 
reflectance),  are  required  to  describe  the  unique  radiance  contribution  of 
the  imaged  water  scene.  However,  in  the  visible  red  region  (and  at  longer 
wavelengths),  water  is  increasingly  opaque  to  the  point  that  only  surface 
reflection  (and/or  emission  in  the  infrared)  requires  description. 
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3)  Except  for  the  rare  instances  of  sustained  windlessness  where 
water  is  smooth  and  specular,  its  surface  is  disturbed  by  the  addition  of 
directional  friction  energy  applied  by  wind  passing  over  it.  Up  to  a  defined 
maximum  wind  velocity,  the  disturbed  surface  remains  analytic  (infinitely 
differentiable)  and  can  be  described  as  a  quasi-stationary,  pseudo- 
Gaussian  process;  the  resulting  power  spectrum  approximations  for 
elevation,  slope,  and  curvature  can  be  analyzed  by  Fourier  methods 
[Kinsman,  1965].  See,  for  example.  Figure  1.3-1. 


Figure  1.3-1.  Schematic  representation  of  the  energy  contained  in  the  surface  waves  of 
the  oceans  -  in  fact,  a  guess  at  the  power  spectrum  [from  Kinsman,  1965]. 


4)  As  an  extension  of  3),  waterwaves  can  be  imaged  over  several 
orders  of  spatial  frequency.  Also,  the  well-behaved  surface  statistics  of 
waterwave  surfaces  (e.g.,  mean  slope  =  0°,  3-sigma  slope  «  30°)  allow 
for  a  broad  range  of  off-vertical  imaging  geometries  without  having  to 
consider  the  complex  effects  of  surface  obscuration.  See,  for  example. 
Figure  1.3-2. 

5)  Also  as  an  extension  of  3),  the  assumption  of  quasi-spatial  and 
quasi-temporal  invariance  of  local  surface  spectra  allows  for  considerable 
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tolerance  in  both  image  collection  and  processing:  precise  pointing  to  an 
exact  ground  reference  point  by  multiple  sensors  is  not  critical;  precise 
collection  timing  by  multiple  sensors  is  not  critical;  and  the  registration  of 
multi-sensor  imagery  is  not  required. 


Figure  1.3-2.  Facet  slope  angle  (beta)  probability  distribution  for  a  wind-roughened  water 
surface,  computed  for  various  wind  speeds  [from  Sidran,  1981], 


6)  Finally,  the  natural  air-water  interface  represents  a  strong 
polarizing  dielectric.  High  contrasts  within  and  between  polarized  images 
should  be  readily  detectable  in  a  significant  range  of  imaging  geometries. 

In  sum,  this  study  will  exploit  one  of  the  more  analytically 
straightforward  terrestrial  surface  features,  and  its  simulation,  for  the 
evaluation  of  polarized  (Stokes  vector)  imagery  in  the  extraction  of  texture 
information  from  more  complex  scenes  of  spatially  correlated  features 
(such  as  cultivated  forests  and  fields)  that  can  be  described  by  their 
spatial  power  spectra. 
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2.0 


Literature  review 


The  full  literature  review  consists  of  three  sections: 

Overview  of  optical  polarization  -  This  section  provides  a 
fundamental  description  of  optical  polarization  and  its 
measurement.  The  work  of  Maxwell,  Stokes,  Malus,  Brewster, 
Fresnel,Young,  Snell,  eta!.,  is  highlighted. 

Polarimetric  image  analysis  -  This  section  provides  a  review  of 
research  that  advances  the  use  of  imaging  techniques  to  measure 
the  polarization  state  of  sensed  spatial  radiation  distributions. 

Spatial  information  content  analysis  -  This  section  provides  a 
review  of  research  that  advances  the  use  of  imaging  techniques  to 
analyze  the  spatial  information  content  of  scenes,  to  include  spatial 
spectrum  analysis.  The  analysis  of  remotely  sensed  waterwave 
scenes  is  emphasized. 

2.1  Overview  of  optical  polarization 

Maxwell's  Equations 

The  property  of  polarization  is  common  to  all  transverse  vector 
waves,  to  include  all  electromagnetic  waves.  The  complete  description  of 
an  electromagnetic  wave  requires  four  field  vectors: 

E  =  the  electric-field  strength, 

B  =  the  magnetic-flux  density. 


6 


D  =  the  electric-displacement  density,  and 

H  =  the  magnetic-field  strength, 


that  are  interrelated  by  Maxwell's  equations: 
D  =  s^so  E , 

B  =  Pmfio  H  , 

div  D  =  47ip, 

div  B  =  0, 

cD/dt  =  W8m8o  dE/dt  =  c  curl  H  -  47tj, 

and 

dB/dt  =  dH/dt  =  -c  curl  E  , 

where 


Equation  2-1 


Equation  2-2 


Equation  2-3 


Equation  2-4 


Equation  2-5 


Equation  2-6 


s  =  the  electric  permittivity,  or  dielectric  constant,  of  the  medium 
under  consideration  (sm)  and  of  free  space  (eo), 
p  =  the  magnetic  permeability  of  the  medium  under 
consideration  (i^m)  and  of  free  space  (i^o). 
p  =  the  electric  charge  density, 
c  =  the  velocity  of  light  in  free  space, 
j  =  the  electric  field  current. 
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w  =  the  cross-sectional  area  of  the  medium  normal  to  the 
direction  of  propagation, 

and  where  the  vector  operators,  div  (divergence)  and  curl,  can  be 
expressed  in  the  form  of  differential  equations: 


Divergence  of  F  =  ^  +  — 

dx  dy  dz 


Equation  2-7 


and 


CURLofF  = 


Equation  2-8 

lydy  d.  ) 

dx  J 

\dx  dy  y J 

Polarization  of  an  electromagnetic  wave  refers  to  the  variation  of 
one  of  these  four  vectors  with  respect  to  time,  t,  at  a  fixed  spatial 
observation  point,  s.  Of  the  four  field  vectors,  convention  defines  the  state 
of  light  polarization  relative  to  the  E-vector;  the  reasoning  behind  this 
choice  of  convention  is  that,  for  light-matter  interactions,  the  electric  field 
force  exerted  on  the  material  electrons  by  a  light  wave  is  much  larger  than 
that  of  its  magnetic  field  force: 


"The  total  force  exerted  by  the  electromagnetic  field  on  a 
particle  of  charge  q  moving  with  velocity  v  consists  of  two 
terms:  the  electric  force  qE  and  the  Lorentz  force  qv  x  B. 

The  ratio  of  magnitudes  of  the  latter  to  the  former  cannot 
exceed  vB/E  or  v/c  where  c  is  the  velocity  of  light.  Because 
v/c  «  1  (for  all  cases  of  interest),  the  Lorentz  force  can  be 
neglected."  [Azzam  &  Bashara,  1977] 

Therefore,  all  subsequent  polarization  analysis  will  limit  the  description  of 
light-matter  interactions  to  the  behavior  of  the  electric-field  strength  vector 
E(s.t). 


8 


Figure  2.1-1.  Curl  relationships 
between  E-  and  H-  vectors  [from 
Partington,  1953]. 


Figure  2.1-2.  E-  and  H-vectors  in  a  linearly 
polarized  monochromatic  wave  [from  Jenkins  & 
White,  1976], 


Stokes  Parameters 


Any  electromagnetic  wave  of  arbitrary  polarization  can  be 
mathematically  represented  by  two  orthogonal  linearly  polarized  complex 
waves; 


E(s,t)  =  E,(z,t)  +  Ey(z,t) . 


Equation  2-9 


assuming  s  specifies  the  coordinates  of  a  right-handed  Cartesian  system 
with  the  wave  propagating  along  the  z-axis  in  the  positive  direction.  The 
two  complex  waves  can  then  be  analyzed  into  their  real  and  imaginary 
components: 


Equation  2-10 

F  p'®  -  F  p'®’'  +  F  p'®y 


=  Exo  [cos(Ox)  -  i  sin(Ox)]  +  Eyo  [cos(Oy)  -  i  sin(Oy)] 

Finally,  the  general  description  of  the  two  real  component  waves  requires 
eight  parameters  as  functions  of  z  &  t: 


Ex(z,t)  =  Exo  cos(Ox)  =  Exo(t)  cos  [(co  t  -  kz)  +  \{i)] 


Equation  2-1 1 


and 

Ey(Z,t)  =  EyO  COS(Oy)  =  Eyo(t)  COS  [( CO  t  -  kz)  +  8y(t)] , 


Equation  2-12 


where 


Ex,Ey  =  the  values  of  the  electric  field  in  the  x  &  y  direction, 
respectively,  at  position  z  and  time  t, 

Exo.Eyo  =  the  amplitude  of  the  electric  field  oscillation  in  the  x  & 
y  direction,  respectively,  at  time  t  and  z  =  0, 

CO  =  the  angular  frequency  of  the  wave, 

k  =  the  wavenumber  of  the  wave  =  27c/wavelength,  and 

6x,5y  =  the  phases  (or  epochs)  of  the  electric  field  oscillation 
in  the  X  &  y  direction,  respectively,  at  time  t. 

Only  four  of  these  parameters  -  Exo,  Eyo,  6x,  and  5y  -  are  required  to 
establish  the  equation  of  the  state  of  the  polarization  ellipse  at  some  fixed 
point  along  the  wavetrain,  e.g.  z  =  0: 

Ex  _L  Ey  ^  Ex  Ey^Q^b  y~b  x)  _  •  2 1 

— 7  +  — 7  -  -  -  sm 

Exo  Ey,  ExoEy, 


Equation  2-1 3 
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Note  that  the  cross-product  exists  because  the  coordinate  axes  do  not,  in 
general,  coincide  with  the  axes  of  the  ellipse. 


Figure  2.1-3.  Specification  of  a  polarization  ellipse,  using  Shurcliffs  notation: 
azimuth  =  a,  ellipticity  =  b/a,  handedness  =  right  (+)  [adapted  from  Shurcliff,  1962]. 


The  Stokes  vector,  containing  the  four  Stokes  parameters,  provides 
a  more  empirically  accessible  description  of  polarization  for  the  more 
general  case  of  partially  polarized  light,  since  they  are  directly  derivable 
from  measurable  irradiance  values: 


S  =  { I,  Q,  U,  V, } , 


Equation  2-14 


where 


I  =  (cs„/2)  [<E,„"(t)  +  E,„^t)>l , 


Equation  2-15 


the  total  irradiance  of  the  light  wave  without  respect  to 
polarization  (always  positive). 


Q  =  (cs„/2)  [<E,„"(t)  -  E,„'(t)>l , 


Equation  2-16 
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the  differential  irradiance  measuring  the  preference  for  linear 
polarization  along  the  x-axis  (positive),  y-axis  (negative),  or 
no  preference  (zero), 


U  =  (cso/2)  [<2  E,o(t)  Eyo(t)  cos(6y(t)  -  5,(t))>] , 


Equation  2-17 


the  differential  irradiance  measuring  the  preference  for  linear 
polarization  along  the  +45°  bisector  (positive),  -45°  bisector 
(negative),  or  no  preference  (zero),  and 


V  =  (C8o/2)  [<2  E,o(t)  Eyo(t)  sin(5y(t)  -  6,(t))>], 


Equation  2-18 


the  differential  irradiance  measuring  the  preference  for 
circular  polarization  that  is  right-handed  (positive),  left- 
handed  (negative),  or  no  preference  (zero). 


Figure  2.1-4  a,b,c,d.  Schematic  representation  of  the  measured  Stokes  vector 
parameters  I,  Q,  U,  and  V,  respectively. 


For  the  special  case  of  monochromatic  light,  which  by  definition  is 
perfectly  polarized,  the  quantities  Exo(t),  Eyo(t),  6x(t),  and  5y(t)  are  time- 
independent,  i.e.,  temporal  averaging  is  not  required. 

<  m  >  =  ^  F(t)  dt 


Equation  2-19 
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Note  that  <  >  specifies  the  temporal  average  of  the  enclosed  quantities. 


The  parameters  of  the  Stokes  vector  provide  a  full  specification  of 
the  polarization  ellipse: 

Scale(Irraaance>^l=(J)[E.;  *  E,:] 


Equation  2-20 


U  2E^E  cos(d 
Azimuth  =  -= 

Q  Ej 


yo 


\v\  sin(6  -6,) 

Shape  =  U  =  ^ 


I  E  ^+E  ^ 

Jto  yo 


and 


Handedness  =  Sign  Of  V 


Equation  2-21 


Equation  2-22 


Equation  2-23 


Another  important  description  is  fully  consistent  with  the  Stokes 
vector  specification; 


Irradiance  =  / 


Equation  2-24 


Degree  of  Polarization  =  P  = 


^Q^  +  U'  +  V^ 
I 


Equation  2-25 


0  —  —  tan’^ 

U 

Equation  2-26 

2 

Q 

and 
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f 


Ellipticity  Angle  ~  s 


f  \ 

1  .  V 

—  sin  /  , 

2  iVe  +t/'+FV 


Equation  2-27 


The  Stokes  vector  can  then  be  recast  using  these  parameters: 


S  =  I  {1,  P  cos(28)  cos(2©),  P  cos(2s)  sin(2©),  P  sin(2£)} 


Equation  2-28 


Alternatively,  a  derivative  Stokes  vector  can  be  specified  using  these  four 
independent  parameters: 


^derivative 


=  {  I,  P,  ©,  s  } 


Equation  2-29 


As  an  introduction  to  the  next  section  of  this  chapter,  it  is  this 
derivative  Stokes  vector  that  is  of  primary  interest  within  this  study.  First, 
these  derivative  Stokes  parameters  have  a  more  direct  correlation  with 
the  physical  properties  of  imaged  features  from  which  reflected  polarized 
radiance  is  measured.  Second,  the  Stokes  parameters  Q  and  U  have  a 
dependent  variation  with  the  azimuth  angle  of  polarization  ellipse  ©,  as 
specified  in  Equation  2-27.  As  will  be  shown  later,  direct  evaluation  of  the 
derivative  ©  (or  T)  image  provides  a  sensitive  indication  of  the  azimuthal 
deviation  in  polarimeter  reference  plane  from  its  expected  alignment  with 
the  principal  plane,  a  deviation  that  would  otherwise  be  masked  as  a 
systematic  error  within  both  the  Q  and  U  images. 
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Figure  2.1-5  provides  examples  of  some  common  polarization  forms  and 
their  Stokes  vector  representations.  Note  that  the  current  Stokes'  notation 
convention  {l,Q,U,V}  is  replaced  here  by  ShurclifTs  notation  {l,M,C,S}, 
after  Perrin  and  Jones.  (Stokes’  initial  notation  was  {A,B,C,D}.) 


Table  2-1 

Stokes  parameter  conventions 


Equation 

Current 

Shurcliff 

Stokes 

[2-14] 

1 

1 

A 

[2-15] 

Q 

M 

B 

[2-16] 

U 

C 

C 

[2-17] 

V 

S 

D 
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Polarization  form  Normalized  Stokes  vector,  Jones  vector 


Sectional  a  y  Standard 

pattern  (deg)  b/a  Ay  I  Ax  (deg)  {/,  M,  C,  51  normalized  Full 


90  0  « 


45  0 


(1,  1,  0,  0} 


(1,  -1,  0,  0} 


{1,  0,  1.  0) 


-45  0  1  ±180  (1,  0,  -1,  0} 


General  linear  0  Any  0 

posi-  or 

tive  ±  180 


{1,  cos  2a!,  sin  2a,  0} 


[  2]  [ 

_[  a  [  2^] 

n  r 
2  L  u  L 

n  r  A,e^-\ 
2  L-iJ 

r  cos  jRi  r  At^^'i 

LlsinjRj 


1,R  1  90  {1,  0,  0,  l! 


—  1,  L  1  -90  {1,  0,  0,  -1} 


0  ~^R  i  90  [1,  0.6,0,  0.81 


1  2rV  5 

90  pi?  2  90  ll,  -0.6,  0,  0.8)  - 

2  5 


r-n  rA,ei‘^  i 

2  L  ij  L^xe‘<‘'+'/2)J 

n  rA,ei^  1 

2  L  ij 


•?ri]  p-"] 

tR]  [::«] 


22.5  0.318,  0.518  45 

R 


General  elliptical 


cos  2(0  cos  2X  r (cos  i2)e'‘*2  |  f  | 

cos  2w  sin  2X  /  ■  • 

sinao  L<=mi?)«‘2  J  U«e«.J 


(flz*  +  flv*) 

General  elliptical,  1  (cf**  —  Uy*) 

partially  polarized  (a,*  +  {2axay  cos  7) 

L{2a,ay  sin  7)^ 


Unpolarized 


{1,  0,  0,  0]  None 


Figure  2.1-5.  Stokes  vector  parameters  for  various  polarization  forms 
[from  Shurcliff,  1962]. 


Reflection  Polarization 

In  1808,  the  French  scientist  Etienne-Louise  Malus  discovered  the 
polarization  of  light  by  its  reflection  from  dielectric  media; 
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"The  Paris  Academy  had  offered  a  prize  for  a  mathematical 
theory  of  double  refraction,  and  Malus  accordingly  undertook 
a  study  of  the  problem.  He  was  standing  at  the  window  of 
his  house  in  the  Rue  d’Enfer  one  evening,  examining  a 
calcite  crystal.  The  Sun  was  setting,  and  its  image  reflected 
toward  him  from  the  windows  of  the  Luxemborg  Palace  not 
far  away.  He  held  up  the  crystal  and  looked  through  it  at  the 
Sun's  reflection.  To  his  astonishment,  he  saw  one  of  the 
double  images  disappear  as  he  rotated  the  calcite.  After  the 
Sun  had  set,  he  continued  to  verify  his  observations  into  the 
night,  using  candlelight  reflected  from  the  surfaces  of  water 
and  glass."  [Hecht,  1987] 

Andrews  [1960]  offers  some  later  details; 


"After  Malus'  discovery  Thomas  Young  pondered  for  eight 
years  how  light  waves  could  be  polarized.  Polarization  was 
being  used  as  evidence  for  the  corpuscular  theory  of  light. 
Suddenly  the  idea  occurred  to  Young  that  if  light  waves  were 
transverse  they  could  be  polarized. 

Sir  David  Brewster  discovered  a  simple  law  to  define 
Malus'  certain  angle  of  incidence  for  which  the  light  is 
reflected  from  the  dielectric  surface  is  totally  plane  polarized: 
The  reflected  light  is  plane  polarized  when  the  angle 
between  the  reflected  ray  and  the  refracted  ray  is  %I2 
[radians  or  90°]." 


This  certain  angle  of  incidence  ©jp  at  a  dielectric  interface  is  now  referred 
to  as  Brewster's  angle.  The  simple  geometric  relation  between  Brewster's 
law  and  Snell's  law  is; 


nj  sin  0jp  =  nt  sin  ©t  (Snell's  Law) 

©p  +  ©t  =  90°  ->  ©t  =  90°  -  ©p 


Equation  2-30 


Equation  2-31 


nj  sin  ©ip  =  nt  sin  (90°  -  ©p)  =  nj  cos  ©p 


Equation  2-32 
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and 


(nt  /  Hj)  =  (sin  0p  /  cos  ©p)  =  tan  0p  (Brewster's  Law) 


Equation  2-33 


Figure  2.1-6  provides 
two  illustrations  of  the 
geometric  relationship 
between  Brewster’s 
law  and  Snell’s  law  at 
a  dielectric  interface. 


Figure  2.1-6  (a)  and  (b).  Two  schematic  representations  of 
the  geometric  relations  between  Brewster’s  law  and  Snell’s 
law  [from  Hecht,  1 987], 
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Application  of  Maxwell's  equations  to  reflection  polarization 


Electromagnetic  theory  can  be  applied  to  confirm  the  empirical  laws 
for  reflection  (to  include  Snell's  and  Brewster's)  and  to  determine  the 
polarization  of  the  reflected  beam,  i.e.,  the  amplitude,  phase,  and  direction 
of  the  electric  field  E. 

The  laws  of  reflection  and  refraction  can  be  derived  from  Maxwell's 
equations  by  imposing  four  boundary  conditions: 

1)  the  tangential  components  of  magnetic  force  H  are  continuous 
across  the  dielectric  boundary  (Ampere's  circuital  law), 

2)  the  tangential  components  of  electric  force  E  are  continuous 
across  the  dielectric  boundary  (Faraday's  law  of  induced 
electromotive  force), 

3)  the  normal  components  of  electric  displacement  D  are 
continuous  across  the  dielectric  boundary  (Gauss'  law),  and 

4)  the  normal  components  of  magnetic  induction  B  are  continuous 
across  the  dielectric  boundary  (Gauss'  law). 

Electromagnetic  theory  predicts  that  when  a  wave  crosses  a  dielectric 
boundary,  i.e.,  a  discontinuity  in  s  and  p,  the  amplitude  reflection 
coefficients  are  dependent  on  the  angle  between  the  electric  vector  and 
the  plane  of  incidence.  These  coefficients  are  determined  from  the 
following  relations: 
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1)  the  vectors  E,  H,  and  the  propagation  vector  P,  are  defined 
through  Poynting's  theorem:  S  =  E  x  H  (x  is  the  cross-product), 


2)  from  the  ratio  E  /  H  =  V(ppo/sso),  refractive  index  n  =  Vs,  p  «  1  for 
both  media  at  optical  wavelengths,  the  relation  between  amplitudes 
E,  H,  and  n  is 


H  =  n  V(so/ppo)  E 


Equation  2-34 


n  Vsq  E  ,  and 


3)  the  components  of  E  and  H  that  are  parallel  to  the  dielectric 
boundary  are  continuous  across  it. 

The  remaining  derivation  follows  from  Clarke  and  Grainger  [1971],  using 
their  figures  and  notation. 


Case  1  ■  E|  rparalleh  to  the  incident  plane  [refer  to  Figure  2.1-7(a)1. 


The  continuous  boundary  condition  for  E  and  H  is  satisfied  by 


Eii  cosx,-  +  Er\\  cosx,  =  E,\\  cosx, 

and 


Equation  2-35 


HiA.  -  HrL  =  Ha 


Equation  2-36 


The  refractive  index  condition  for  E  and  H  is  satisfied  by 


wEii  -  niEr\\  -  n2Ei\ 


Equation  2-37 


"The  amplitude  reflection  coefficient  is  Er||  /  Ej||  =  ry,  say,  and  is 
easily  seen  from  [the  equations  above]  and  the  law  of  reflection  (X|  =  Xr),  to 
be  given  by: 


n2 cosx-  -  mcosx, 

w^cosx,  +  mcosx, 


Equation  2-38 


Using  the  law  of  refraction  that  n^  sin  Xj  =  r\2  sin  Xj,  [this  equation] 
reduces  to: 

tan(^X,  -  X,) 
tan(Xi  +  X,) 


Equation  2-39 


For  angles  of  incidence  not  very  different  from  zero,  the  two 
tangents  are  both  positive  and  so  the  minus  sign  indicates  that,  in  fact,  the 
reflected  wave  undergoes  a  phase  shift  of  n  [radians]  with  respect  to  the 
situation  shown  in  [figure  (a)],  i.e.  nodality  is  predicted.  For  (Xj  +  Xt)  =  nl2 
[radians],  the  reflection  coefficient  for  Ey  falls  to  zero  and  so  the  reflected 
beam  can  only  contain  a  perpendicularly  vibrating  component.  It  is  easy  to 
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show  that  (Xj  +  xO  =  7i/2  [radians]  corresponds  to  a  value  for  Xj  of 
arctan(n2/ni).  This  angle  of  incidence  is  known  as  Brewster's  angle.  For  X| 
greater  than  Brewster's  angle,  r||  changes  sign  and  becomes  positive." 
[Clarke  &  Grainger,  1971] 

Case  2.  Ex  (perpendicular^  to  the  incident  plane  [refer  to  Figure  2.1-7(b1]. 


The  continuous  boundary  condition  for  E  and  H  is  satisfied  by 


Hii  cosx,-  +  Hri  cosx,  =  Hii  cosx, 


Equation  2-40 


and 


EiX  -  Erl  -  Eil 


Equation  2-41 


The  refractive  index  condition  for  E  and  H  is  satisfied  by 


niEiiCosXi  -  niEriCOSX^  =  m^acosx, 


Equation  2-42 


The  amplitude  reflection  coefficient  is  Er±  /  Ejx  =  rx;  once  again, 
combining  the  equations  above  with  the  law  of  reflection  (X|  =  Xr),  the 
reflection  coefficient  for  Ex  is: 


=  mcosx,  -  mcosx, 
mcosx,  +  n2^osx, 

Again,  using  the  law  of  refraction  that  n^  sin  X|  =  n2  sin  Xj,  this 
equation  reduces  to: 


Equation  2-43 
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Equation  2-44 

sin(Xi  +  X,)  I - - 

The  key  insight  in  reflection  polarization  is  that,  at  Brewster's  angle, 
energy  could  only  be  directed  to  the  production  of  an  incident-plane  wave 
if  the  wave  was  a  longitudinal  wave.  The  observation  by  Brewster, 

Fresnel,  and  ultimately  Young  (in  1817)  that  no  such  polarized  wave  is 
produced  was  taken  as  proof  that  light  vibrations  are  transverse. 


The  reflection  matrix  defined  by  Mueller  calculus  is  of  the  form: 


Mueller  calculus  will  predict  the  Stokes  vector  of  the  reflected  light  by  pre¬ 
multiplying  the  reflection  Mueller  matrix  R  by  the  Stokes  vector  of  the 
incident  light  ray  S: 

^reflection  “  ^  ^incident  Equation  2-46 

The  Mueller  calculus,  based  on  the  work  of  Hans  Mueller,  is  an  empirical 
matrix-based  formalism  for  treating  all  the  possible  physical  transforma¬ 
tions  of  Stokes  vectors.  Since  the  Stokes  vector  is  a  four-parameter 
vector,  the  Mueller  operators  are  all  4  x  4  matrices. 
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Without  resorting  to  the  inclusion  of  a  coordinate  transformation  matrix, 
the  Stokes  parameters  here  are  presumed  to  be  measured  such  that  the 
z-axis  corresponds  to  the  direction  of  propagation  and  the  x-  and  y-axes 
are  respectively  perpendicular  and  parallel  to  the  plane  of  incidence. 

The  benefits  of  using  this  type  of  formalism  (the  Mueller  calculus)  are 
noted  by  Clarke  and  Grainger  [1971]: 


"This  one  matrix  succinctly  summarizes  all  the  behaviour  we 
have  been  describing,  including  the  distinction  between 
phase  changes  and  nodality.  For  example,  the  [V]  Stokes 
parameter  of  the  reflected  light  will  have  the  opposite  sign  to 
that  of  the  incident  light  if  r±  and  rj  have  the  same  sign,  as 
they  do  for  angles  of  incidence  less  than  Brewster's  angle. 
Thus  the  matrix  predicts  a  change  of  handedness  under 
these  circumstances.  For  angles  of  incidence  greater  than 
Brewster's  angle,  rx  and  rj  have  opposite  signs,  and  hence 
the  [V]  parameters  of  both  incident  and  reflected  light  have 
the  same  sign  -  indicating  no  handedness  reversal. 

Another  interesting  prediction  of  this  matrix  is  that  the 
azimuths  of  the  incident  and  reflected  beams  are,  in  general, 
different.  This  is  well  illustrated  by  considering  the  reflection 
of  linearly  polarized  light  at  normal  incidence,  when  rx  =  rj  = 
r.  The  matrix  becomes: 

'10  0  O' 

,010  0 
R  =  r^ 

0  0-10 

0  0  0  -1 


Application  of  this  matrix  to  the  Stokes  vector,  {[l,Q,U,V]},  of 
the  incident  light,  shows  that  the  reflected  light  is  linearly 
polarized,  but  with  the  opposite  sign  of  [U  /  Q].  The  azimuths 
of  the  incident  and  reflected  light  are  thus  not  the  same." 
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2.2  Optical  polarimetric  imaging 


Overview 

Optical  imaging  polarimetry  is  a  relatively  new  measurement 
technique.  Although  the  required  methods  and  materials  have  been 
available  for  some  time,  the  first  demonstrations  and  applications  of 
optical  imaging  polarimetry  did  not  appear  in  literature  until  the  1970s. 

This  unique  research  activity  is  distinguished  from  the  substantial  quantity 
of  current  studies  that  utilize  non-imaging  optical  polarimetric  techniques, 
e.g.,  photopolarimetry:  and  also  the  substantial  number  of  current  studies 
that  utilize  non-optical  imaging  polarimetric  techniques,  e.g.,  polarimetric 
synthetic  aperature  radar  (SAR). 

Initial  studies  evaluated  the  feasibility  of  optical  imaging  polarimetry 
in  the  enhanced  discrimination  and  classification  of  terrestrial  features,  in 
particular,  vegetation  [Curran,  1981b,  1982;  Walraven,  1981;  Egan,  1985; 
Duggin  etal.,  1989;  Egan  etal.,  1992]  and  soil  moisture  [Curran,  1978, 
1979,  1980,  1981a;  Egan,  1985].  These  studies  utilized  a  single  35-mm 
camera  with  a  linear  polarizer  that  was  manually  rotated  90°  in  order  to 
achieve  discrimination  of  the  first  two  Stokes  parameters  and  obtain  an 
estimate  of  the  degree  of  polarization  (Q/l). 

Prosch  etal.  [1983]  constructed  a  video  imaging  polarimeter to 
study  the  polarization  of  solar  radiation  reflected  from  the  natural 
environment.  The  system  was  composed  of  three  vidicons  aligned  in 
parallel,  with  linear  polarizers  rotated  at  0°,  60°,  and  120°  relative  to  the 
reference  plane  such  that  the  Stokes  parameters  I,  Q,  and  U  could  be 
measured.  The  composite  color  video  signal  was  coded  so  that  deviation 
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from  white  would  represent  increasing  polarization;  tint  specified  the 
azimuth  angle  of  the  polarization  ellipse  whereas  hue  defined  the  degree 
of  polarization. 

Designs  for  a  spaceborne  imaging  multispectral  polarimeter  have 
been  proposed  [Egan,  1986].  While  a  full  scale  testbed  was  infeasible, 
limited  spaceborne  polarimetric  photography  experiments  were  executed 
onboard  six  Space  Shuttle  missions,  STS  51 A  (October  1984)  through 
STS  28  (August  1989)  [Coulson  etal,  1986;  Duggin  eta!.,  1989;  Egan  et 
al.,  1991,  1992].  In  these  experiments,  two  70-mm  Hasselblad  cameras 
were  aligned  in  parallel,  with  linear  polarizers  (Polaroid  HN22)  rotated  at 
0°  and  90°  relative  to  a  reference  plane:  only  the  Stokes  parameters  I  and 
Q  were  measured.  A  significant  consideration  for  accurate  measurement 
was  the  ability  of  the  mission  specialist  to  align  the  two-camera  system 
with  the  mean  principal  plane  of  the  imaged  radiance  field.  Also,  the 
experiment  suffered  from  the  lack  of  sensitometric  calibration,  particularly 
important  for  a  dual  camera  (hence  dual  film  roll)  system.  Again,  the  first 
two  Stokes  parameters  provided  an  estimate  of  the  degree  of  polarization 
as  a  simple  ratio  (Q/l). 

Based  upon  the  promising  results  of  the  Space  Shuttle  experiment, 
the  US  Army  flew  a  helicopter-borne  video  polarimeter  over  various  cover 
types  and  spanning  a  range  of  imaging  geometries  [Israel,  1991].  In  this 
experiment,  a  rotating  linear  polarizer  was  mounted  on  a  single  video 
camera:  each  imaged  scene  collection  contained  a  complete  360°  rotation 
of  the  polarizer. 
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Figure  2.2-1  illustrates  the 
execution  of  the  Space 
Shuttle  polarimeter  exper¬ 
iment  onboard  STS  51A. 
One  of  the  concerns  in  the 
evaluation  of  results  was 
the  extent  of  depolariza¬ 
tion  created  by  imaging 
through  the  multi-layer 
laminate  shuttle  windows 
[Egan  et  a/., 1991]. 


Figure  2.2-1.  NASA-JSC  photograph  5 1A-S 19-08-008 
[from  Israel,  1991). 


The  most  current  operational  imaging  polarimeter  is  the  airborne 
POLDER  (Polarization  and  Directionality  of  the  Earth’s  Reflectances) 
instrument,  first  flown  in  1990  under  the  auspices  of  the  French  space 
agency,  I’Centre  National  d’Etudes  Spatiales  [Breon  &  Deschamps,  1993; 
Deuze  etal.,  1993;  Deschamps  etai,  1994].  The  stated  mission  of 
POLDER  is  to  measure  the  spectral  bi-directional  reflectance  and 
polarization  of  solar  radiation  reflected  by  the  earth-atmosphere  system. 
The  instrument  consists  of  a  CCD  detector  array,  a  rotating  wheel 
containing  16  filter  slots,  and  a  wide-field-of-view  telecentric  lens  system. 
One  slot  on  the  rotating  wheel  contains  an  opaque  filter  for  measuring 
detector  dark  current,  six  slots  carry  unpolarized  spectral  filters,  and  nine 
slots  carry  polarized  spectral  filters  (three  spectral  regions  with  0°,  60°, 
and  120°-rotated  linear  polarizations  for  each  spectral  region). 

Figure  2.2-2  illustrates  the  main  components  of  the  POLDER  instrument. 
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1994]. 

Deschamps  et  al.  [1994]  describe  plans  to  place  a  modified 
POLDER  sensor  onboard  the  Japanese  ADEOS  (Advanced  Earth 
Observing  Satellite),  scheduled  for  launch  in  1996,  in  support  of  the  World 
Climate  Research  Program  (WCRP)  and  the  International  Geosphere  and 
Biosphere  Program  (IGBP).  The  only  other  planned  spaceborne  system 
capable  of  measuring  global  polarized  reflectance  is  the  Earth  Observing 
Scanning  Polarimeter  (EOSP),  scheduled  for  launch  in  2003  [Travis, 
1992], 


In  sum,  optical  polarimetric  imaging  for  remote  sensing  purposes  is 
a  recently  emergent  technology,  with  a  modest  level  of  investigation  to- 
date  compared  with  the  polarimetric  studies  of  active  imaging  systems, 
e.g.,  SAR  and  LASER. 

Review  of  Walraven  [1977.  1981] 

A  detailed  review  of  Walraven’s  [1977,  1981]  polarimetric  imaging 
technique  follows.  The  main  purpose  of  this  review  is  to  present  the 
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underlying  Stokes  parameter  derivation,  since  the  current  study  duplicates 
much  of  Walraven’s  basic  approach  to  Stokes  parameter  measurement. 

An  important  element  of  Walraven’s  and  preceeding  remote 
sensing  investigations  is  the  deliberate  omission  of  the  measurement  of 
the  Stokes  parameter  V,  the  magnitude  of  circular  polarization.  For 
terrestrial  imaging,  the  V  Stokes  parameter  of  the  reflected  light  is  two  to 
three  orders  of  magnitude  smaller  than  Q  and  U,  is  therefore  assumed  to 
be  near-zero,  and  absorbed  within  the  overall  system  error  [Talmage  & 
Curran,  1986].  The  simplifying  benefit  of  this  omission  is  that  polarimeter 
design  and  construction  requires  the  use  of  linear  polarizers  only. 

The  generalized  Mueller  matrix  description  for  an  ideal  linear 
polarizer  is  of  the  form: 
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sin  20 


cos  20 
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cos  20  sin  20 


sin  20 
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Equation  2-47 


where  0  is  the  angle  of  rotation  of  the  polarizer  with  respect  to  the 
reference  axis.  The  resulting  transformation  of  an  incident  Stokes  vector 
by  an  ideal  linear  polarizer  can  then  be  described  by  Mueller  calculus: 


=  M,S 


incident 


Equation  2-48 


Again,  the  Mueller  calculus  is  an  empirical  matrix-based  formalism  for 
treating  the  physical  transformation  of  radiation  as  described  by  Stokes 
vectors;  in  this  case,  it  treats  the  effect  of  transmission  through  a  linear 
polarizer  through  premultiplication  by  the  matrix  operator  M. 
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Walraven  used  a  single  35-mm  camera  with  a  linear  polarizer.  He 
exposed  four  images  of  each  scene,  rotating  the  polarizer  for  each  image 
at  0°,  45°,  90°,  and  135°,  respectively.  The  effective  Stokes  vector 
transformation  for  each  of  the  image  collections  follows  directly  from 
Equations  2-47  and  2-48  for  0  values  of  0°,  45°,  90°,  and  135°  applied  as 
arguments  to  the  general  Mueller  matrix  : 
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Equation  2-52 


Figure  2.2-3  illustrates  the  preferred  linear  polarization  state  of  the 
transmitted  light  for  each  of  the  four  polarizers. 
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Figure  2.2-3.  Schematic  representation  of  the  four  imaging  measurements  of  the  linear 
polarization  state  of  the  incident  radiance,  after  Walraven  [1981]. 


However,  since  the  film  records  only  the  intensity  I  of  the  filtered  radiance, 
this  transformation  can  be  simplified  to  consider  only  the  effect  of  polarizer 
rotation  on  the  first  Stokes  parameter: 


h=\  [4  +  fe  cos  20)+  (u„  sin  20)] 

For  each  of  these  cases,  the  intensity  of  the  transformed  Stokes  vector 
that  is  measured  on  film  is: 


Equation  2-53 


^  =  ^[4  +  8,.] 
4  =^[4  +  ;^.] 
4  ~  2  \f'"  ~  ] 

and 

45=^[4-4„] 


Equation  2-54 


Equation  2-55 
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Equation  2-57 
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The  recovery  of  the  first  three  parameters  of  the  incident  Stokes 
vector  image  Sjn  is  accomplished  through  addition  and  subtraction  of  the 
four  filtered  images: 


hn  ~  [^0  -^90  ]  “  [^45  A35  ]  ~  ^  [^0  hi  "*■  ^90  hn  ] 


Equation  2-58 


Qin  ~  -^O  ] 


Equation  2-59 
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Equation  2-60 


Equation  2-61 


Consistent  with  Walraven’s  approach,  the  current  study  employs  a 
derivative  transformation  of  the  Stokes  vector  that  has  a  more  direct 
correlation  with  the  physical  properties  of  the  imaged  features  from  which 
the  polarized  radiance  is  reflected: 


S  in  = 
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where 


Equation  2-62 
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4=  the  irradiance  of  the  incident  radiance  (unchanged), 

4=  the  degree  of  polarization  of  the  incident  radiance, 

0;„  =  the  azimuth  angle  of  the  polarization  ellipse  of  the  incident  radiance, 
and 

S;„=  the  ellipticity  angle  of  the  polarization  ellipse  of  the  incident  radiance. 

The  definition  of  these  parameters  is  reprised  here  from  Equations  2-25 
through  2-27: 


+U^  +V^'(^  +U^ 
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Equation  2-63 
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Equation  2-64 
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Equation  2-65 

For  linearly  polarized  light,  the  polarization  ellipse  approaches  a  line 
segment,  with  ellipticity  (b/a)  and  ellipticity  angle  z  both  approaching 
zero. 


The  four  photographic  slides  from  each  imaging  collection  were 
then  digitized  into  a  512  x  512  pixel  array,  using  a  vidicon  camera  with  8- 
bit  (0-255  grayscale)  resolution.  The  four  images  were  then  co-registered 
and  the  derivative  Stokes  parameter  images  I,  P,  and  0  were  calculated 
from  the  application  of  Equations  2-58  through  2-64. 
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A  significant  enhancement  of  the  current  study  over  Walraven’s 
technique  is  the  incorporation  of  sensitometric  co-calibration  prior  to  the 
calculation  of  the  Stokes  parameters.  The  co-calibration,  or  joint 
calibration  of  the  four  images  to  a  single  common  characteristic  curve, 
corrects  for  the  non-linear  transformation  of  light  irradiance  to  emulsion 
density  that  is  inherent  to  film,  especially  in  the  highest  and  lowest 
exposure  regions.  The  effect  of  this  co-calibration  is  to  provide  a  more 
accurate  estimate  of  the  derived  Stokes  parameters  that  are  directly 
calculated  from  measured  irradiance. 
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2.3  Spatial  information  content  analysis  of  imagery 


Introduction 


Within  this  study,  spatial  information  content  analysis  refers  to 
those  analytic  techniques  that  are  applied  to  the  measurement  of  spatial 
scene  structure  -  in  particular,  the  measurement  of  remotely  sensed 
surface  structure  -  through  the  analysis  of  their  spatial  reflected  intensity 
distributions.  This  unique  research  activity  is  distinguished  from  existing 
studies  that  use  spatial  information  content  analysis  to  characterize  the 
performance  of  imaging  systems  independent  of  image  content,  e.g., 
optical  transfer  function  (OTF)  measurement;  and  also  existing  studies 
that  use  spatial  information  content  analysis  to  characterize  surface 
structure  under  controlled  conditions,  e.g.,  LASER  scatterometry. 


Terrestrial  remote  sensing  of  waterwave  surfaces 

In  1925,  Schumacher  made  simultaneous  near-horizontal  stereo¬ 
pairs  from  a  ship  with  the  intent  of  measuring  the  variability  of  wave 
heights.  The  utility  of  this  method  was  severely  limited  due  to  many 
factors:  1)  the  camera  baseline  was  restricted  to  the  length  of  the  ship 
(one  camera  fore  and  aft);  2)  waves  in  the  foreground  obstructed  waves  in 
the  background;  3)  backsides  of  waves  were  not  visible;  and  4)  there  was 
a  lack  of  'ground'  control  on  the  open  seas  for  height  determination  -  the 
errors  are  especially  pronounced  from  an  oblique  perspective  [Pos,1988]. 
This  experiment  represents  a  limiting  case  for  spatial  feature 
measurement  via  image  analysis. 
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In  1933,  Hulburt  [1934]  made  polarized  and  unpolarized  oblique 
photographs  of  sun  glitter  on  sea  waves  to  measure  the  polarization  of 
light  at  sea  with  respect  to  surface  roughness,  sun  angle,  and  weather 
conditions.  Because  the  widths  of  glitter  patterns  correlate  to  the 
maximum  slope  of  the  sea  surface,  Hulburt  was  able  to  demonstrate  that 
waves  in  the  North  Atlantic  varied  from  15°  inclination  when  winds  were 
blowing  at  3  knots  up  to  25°  inclination  at  18  knots. 

Sawyer  [1949]  mounted  a  Sonne  strip  camera  on  a  fast  low-flying 
airplane  to  photograph  narrow  strips  of  the  sea  surface  and  measure  the 
directional  spectrum  of  the  waves.  The  field  of  view  was  too  narrow  to 
capture  significant  amounts  of  surface  data  orthogonal  to  the  flightline. 
The  accuracy  of  the  spectral  estimate  decreased  with  the  angle  from  the 
flightline. 

Also  in  1949,  Barber  [1949,  1954]  analyzed  single  photographs  of 
sea  surfaces  to  determine  wave  direction,  but  he  was  unable  to  determine 
the  two-dimensional  spatial  spectrum  due  to  the  computing  limitations  of 
the  time. 

At  this  point,  it  became  apparent  that  an  essential  requirement  for 
future  analysis  of  spatially  correlated  surfaces  is  to  have  near-nadir,  high- 
resolution,  high-contrast  images  covering  large  areas.  The  intent  is  to 
photographically  capture  significant  information  about  the  largest  range  of 
spatial  frequency  components  without  perspective  distortion  or  hidden 
surface  detail.  It  also  became  apparent  that  the  analysis  of  large  amounts 
of  spatial  information  would  require  the  forthcoming  capabilities  of  the 
digital  computer. 
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In  1951,  Cox  and  Munk  [1954a,  b]  measured  the  wave-slope 
distributions  of  the  sea  surface  from  aerial  photographs  of  sun  glitter 
patterns.  They  computed  the  distribution  from  the  measured  variation  of 
radiance  within  a  glitter  pattern  instead  of  computing  maxima  from  the 
pattern  boundaries  as  done  by  Hulburt.  Four  cameras  were  flown  from  a 
single  airplane  at  altitudes  of  2000  feet,  with  two  used  as  imagers  and  two 
used  as  radiometers.  Their  image  analysis  was  quite  sophisticated;  it 
accounted  for  sun  diameter,  angular  reflectivity,  lens  falloff  (vignetting), 
film  sensitivity,  exposure  calibration  (sensitometry),  and  ultimately 
provided  a  first-order  relationship  between  film  density  and  directional 
wave-slope  probability. 

In  1953,  Schooley  [1954]  performed  a  simplified  version  of  the  Cox 
&  Munk  experiment  by  taking  flash  photographs  of  a  river  surface  from  a 
45-foot  bridge  elevation  at  night.  The  main  limitation  of  his  experiment 
was  the  probable  inhomogeneity  of  the  water  surface  due  to  limited  fetch 
(the  surface  area  where  waterwaves  are  being  generated  by  the  wind) 
and  the  presence  of  wave-refracting  obstacles  in  the  water. 

In  1954,  medium-altitude  (3000  feet)  stereophotography  was 
employed  by  Marks  and  Ronne  [1955]  to  generate  stereopairs  of  sea 
surfaces.  Two  airplanes  carried  radio-synchronized  cameras  and  a 
surface  ship  acted  as  'ground'  control  in  the  photographs.  Elevations  were 
photogrammetrically  measured  at  discrete  points  and  the  sampled 
elevation  array  was  then  autocorrelated  (the  sampling  distance 
determined  the  desired  spatial  resolution).  This  experiment  marks  the  first 
recorded  use  of  a  digital  computer  to  calculate  the  two-dimensional 
spectra  of  waterwaves.  The  work  of  Cote  et  al.  [1960]  enhanced  this  basic 
technique.  More  recent  stereophotogrammetric  efforts  include  Holthuijsen 
[1983a,  b]  and  Pos  etal.  [1988].  Elements  of  this  later  work  include 
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methods  to  render  the  water  opaque  so  that  a  more  exact  calculation  of 
the  height  field  can  be  made. 

During  the  1950s  and  early  1960s,  Longuet-Higgins  [1952-1962; 
Cartwright  &  Longuet-Higggins,1956]  elaborated  on  the  results  of  Cox  and 
Munk  to  formulate  the  statistical  theory  of  patterns,  paths,  number, 
frequency,  and  distributions  of  specular  reflection  points  on  randomly 
moving  surfaces.  Stilwell  [1969;  Stilwell  &  Pilon,  1974]  correlated  the 
statistics  of  sea-surface  images  to  the  wave-slope  statistics  of  the  actual 
sea  surface.  Under  the  assumptions  of  uniform  sky  radiance,  optimized 
viewing  geometry,  and  small  surface  slopes,  Stilwell  derived  the 
relationship  between  the  film  transmittance  of  an  imaged  surface  point 
and  the  range  component  of  the  wave  slope  at  that  surface  point.  He 
further  demonstrated  a  ‘linearizable’  relationship  between  the  spatial 
image  spectrum  and  the  imaged  surface  slope  spectrum.  Kasevitch  [1975] 
extended  Stilwell's  model  to  second  order  to  develop  an  optimization 
criterion  for  the  relationship.  Chapman  and  Irani  [1981]  took  this  work  one 
step  further  by  creating  a  synthetic  model  and  executing  a  limited 
quantification  of  the  error  magnitudes  associated  with  the  parametric 
dependence  of  this  linear  model.  North  [1989]  enhanced  Chapman  and 
Irani’s  model  to  consider  sub-resolution  wave  slopes  and  then  executed  a 
more  comprehensive  parametric  surface  exploration. 

Shores  [1980,  1981]  developed  a  novel  technique  for  remotely 
sensing  surface-flow  velocities  based  on  imagery  of  monochromatic 
wavetrains  of  known  frequency  (such  as  those  generated  by  a  motor  boat) 
propagating  over  the  region  of  interest.  His  work  demonstrated  that  the 
wavelength  and  direction  of  two  wavetrains  provided  all  the  information 
required  to  calculate  surface  flows.  Gotwols  and  Irani  [1980]  developed  a 
similar  technique  to  determine  the  phase  velocity  of  short  gravity  waves. 
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Exotic  sensors  using  LASER  [Palm  etal.,  1977;  Schau,  1978; 
Abshire  &  McGarry,  1987]  and  LIDAR  [Weinman,  1988]  have  been  used  to 
extract  directional  spectra  and  surface  backscatter  data  at  higher  wave- 
numbers  (i.e.,  the  capillary  wave  regime).  Synthetic  Aperture  Radar  (SAR) 
imagery  has  been  used  to  estimate  spectra,  phase  velocities,  and 
propagation  directions  at  lower  wavenumbers  (i.e.,  the  gravity  wave 
regime)  [Monaldo  &  Lyzenga,  1986;  Monaldo  &  Kasevich,  1982;  Carlson, 
1984].  Also,  Long  Wave  Infrared  (LWIR)  sensors  have  been  used  to 
calculate  spatial  spectra  of  ocean-surface  temperature  [Saunders,  1967, 
1968;  McLeish,  1970]. 

Lybanon  [1985]  reported  on  the  implementation  of  an  automated 
image-analysis  system  by  the  U.S.  Naval  Ocean  Research  and 
Development  Activity  (NORDA,  now  the  Naval  Research  Laboratory  or 
NRL).  The  Interactive  Digital  Satellite  Image  Processing  System  (IDSIPS) 
can  automatically  derive  the  sea-surface  slope  statistics  from  sun  glitter 
images  through  analysis  of  the  imaging  geometry.  As  a  late  practical 
example,  Fisher  [1986]  analyzed  four  sun  glitter  images  taken  from  the 
space  shuttle  Challenger  (STS-41G)  to  locate  acoustically  important 
oceanographic  features  in  support  of  hydro-acoustical  sensor  placement. 

Breon  and  Deschamps  [1993]  described  the  use  of  the  POLDER 
(Polarization  and  Directionality  of  the  Earth  Reflectance)  instrument  to 
derive  ocean  wave  slope  distributions  from  polarized  directional  specular 
reflectance  measurements.  Inversion  of  an  analytical  model  against  their 
data  accurately  fit  (<  1 .2%)  all  spectral  and  angular  reflectance  variations 
observed  with  the  POLDER  sensor.  These  observations  included  a  large 
reflectance  asymmetry  relative  to  the  principal  plane  that  they  correlated 
with  anisotropic  wave  slope  distribution  created  by  directional  wind  stress. 
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Spatial  spectrum  analysis  of  correlated  scenes 


A  comprehensive  review  would  include  the  work  of  Stillwell  and 
Pilon  [1974],  and  Kasevich  etal.  [1971,  1972],  The  emphasis  of  this 
earlier  work  is  on  the  analysis  of  coherent  optical  processing  techniques 
as  applied  to  photographic  emulsions  of  spatially  correlated  scenes. 
Kasevich  [1975]  provides  a  general  introduction  to  the  first-order  theory 
subsequent  to  the  development  of  an  geometric-optics-approximation 
second-order  theory  that  estimates  the  optimum  viewing  geometry  for 
obtaining  reasonable  spectra.  Stillwell  [1969]  provides  additional 
development  of  optical  imaging  theory  subsequent  to  performing  an 
optical  analysis  to  derive  two-dimensional  spectra. 

Review  of  Kasevich  [1975] 

The  essential  requirement  for  the  determination  of  spectra  from 
spatially  correlated  scenes  (in  this  case,  waterwave  scenes)  is  to  have  the 
spatial  modulation  of  the  imaged  scene  be  proportional  to  the  surface 
profile.  Kasevich  uses  the  example  of  a  transparency  with  film  exposure, 
E,  defined  over  its  linear  region  by 


E{y)  =  Uyr'^ 


ri+/(T)i 

-Y/2 

/o(t)  . 

Equation  2-66 

where 


E{y)  =  fo{y)  +  fiy) 

and 


Equation  2-67 
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Y  =  the  film  gamma 

such  that 

foiy)  -  the  mean  exposure  on  film 

and 

/ (y)  =  the  exposure  modulation  due  to  scattering  of  radiance  from 
specular  wave-slope  facets. 

This  example  is  given  for  the  one-dimensional  case. 

If  fo(y)  f(y).  fhsn  Equation  2-66  can  be  expanded  in  a  binomial 

series  to  yield  the  approximation 


E(y)  =  fo(yr'" 

f(y)  ^ 

Equation  2-68 

\2J 

\fo(y)J 

The  estimation  of  the  slope  spectrum  from  the  forward  Fourier  transform 
of  the  image  requires  that  f(y)  be  linear  with  respect  to  the  wave  slope 
dz/dy,  where  z  is  the  surface  elevation.  This  condition  can  only  be 
approximately  satisfied  for  waterwave  surfaces  because  of  the  non¬ 
linearity  of  1)  the  spatial  radiance  distributions  found  in  nature,  2)  the 
Fresnel  reflectivity  variation  with  respect  to  incidence  angle,  and  3)  the 
refraction  of  upwelling  subsurface  radiance  in  the  direction  of  the 
observer. 
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Review  of  Stillwell  [1969] 


For  a  small  wave-slope  angle  p  ,  the  small  angle  approximation  is: 


P  -  tan'^ 


dz 

dy^ 


dz 

dy 


Equation  2-69 


where  p  is  the  fundamental  parameter  for  extracting  wave-slope  spectra 
from  imagery.  The  reflected  radiance  observed  at  azimuth  angleG  is  a 
simple  function  of  Fresnel  reflectivity  and  incident  radiance: 


=  L(\x)  R((f)) 


Equation  2-70 


where 


Zo(6,  P,co )  =  the  observed  reflected  radiance, 

Z(p)  =  the  incident  radiance  to  be  reflected, 

R((o )  =  the  Fresnel  reflectivity  (for  any  arbitrary  polarization), 
0  =  the  zenith  angle  of  observation, 

P  =  the  slope  of  the  reflecting  surface  facet, 
p  =  the  zenith  angle  of  the  incident  radiance, 

and 
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CO 


=  the  angle  of  incidence, 


such  that 


(0  =0  -  (3 


and 


|a,  =(0  -  (3 


Equation  2-71 


Equation  2-72 


Figure  2.3-1  illustrates  the  angular  relations  for  the  one-dimensional  case. 


Figure  2.3-1.  One-dimensional  angle  relationships  [from  North,  1989]. 


The  variation  of  observed  radiance  with  respect  to  a  change  of 
surface  slope  c;f3  at  some  point  is 


c/Z(9) 


dL{\x,) 

(i?(co)) 

d\i 

+  (I(p)) 

'dR{&) 

d(£) 

Equation  2-73 

d\i 

_  d(o  _ 

If  (3  is  a  small  angle  and®  s  |j,  ,  then 


=  /■iY^  j  Sl'ci)  J  +  L(t^)R'(a)J^. 


where  the  prime  (')  denotes  the  first  derivative  with  respect  to  the 
argument. 

With  the  assumption  that  the  water  surface  remains  analytic  (i.e.,  infinitely 
differentiable),  the  small-angle  linear  approximation  holds  in  the  Fourier 
transform  for  waveslope  angles  (p )  up  to  at  least  30°. 


Equation  2-74 
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3.0  Methods  and  Materials 


The  full  description  of  methods  and  materials  comprises  four 
sections,  one  for  each  of  the  four  phases  of  this  study: 

Image  Data  Collection  -  This  section  describes  the  design  and 
construction  of  an  imaging  polarimeter  system  and  the  execution  of 
both  in-field  and  aerial  image  data  collections. 

Image  Data  Processing  -  This  section  describes  the  post-collection 
processing  of  the  image  data:  film  development,  image  digitization, 
digital  image  data  manipulation,  and  production  of  Stokes  imagery. 

Image  Data  Simulation  -  This  section  describes  the  model  that  was 
used  to  simulate  Stokes  imagery  of  the  reflected,  polarized 
skydome  under  clear  sky  conditions  -  as  a  basis  for  comparison 
with  an  image  subset  from  the  in-field  data  collection. 

Image  Data  Analysis  -  This  section  describes  the  analytic 
techniques  that  were  applied  to  the  Stokes  imagery,  to  include  two- 
dimensional  fast  Fourier  transformation  (FFT)  of  selected 
subimages. 

In  addition,  the  source  code  developed  in  this  study  for  the  image  data 
simulation  is  found  in  Appendix  A. 
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3.1  Image  Data  Collection 


3.1.1  General 

Construction  of  the  imaging  polarimeter 

A  filter  mount  was  constructed  to  fit  over  the  30-mm  four-lens 
system  of  a  Nishika  N8000  3-D  camera.  Four  small  (18.5  mm  x  18.5  mm) 
square  sections  were  cut  from  the  same  sheet  of  Polaroid  HN38  linear 
polarizer  such  that  the  transmission  axes  of  first  filter  pair  are  rotated  0° 
and  90°  about  the  optic  axis  with  respect  to  a  common  reference  plane; 
and  that  the  transmission  axes  of  the  second  filter  pair  are  rotated  45°  and 
1 35°  about  the  optic  axis  with  respect  to  the  same  reference  plane.  The 
four  filter  squares  were  aligned  in  the  mount,  secured  in  place,  and  the  full 
assembly  was  then  fitted  over  the  four-lens  system  and  secured.  If  the 
four  small  filter  squares  were  cut  perfectly,  the  abutment  of  their  edges 
against  each  other  and  against  the  interior  edge  of  the  mount  would 
ensure  rotational  alignment  of  the  linear 
polarizers. 

Figure  3.1-1  illustrates  the  resulting 
imaging  polarimeter  system.  The  mount 
containing  the  four  rotated  linear  filters  lies 
between  the  camera  and  two  samples  of 
Polaroid  HN38  linear  sheet  polarizer. 

Alignment  of  the  imaging  polarimeter 

A  key  assumption  is  that  the  optic 
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axes  in  the  four-lens  system  are  all  in  parallel;  by  definition,  the  common 
reference  plane  for  filter  rotation  contains  these  four  optic  axes.  Another 
key  assumption  is  that  the  transmission  axis  of  the  HN38  polarizer  is  in 
parallel  with  the  specified  edge  of  the  manufacturer-supplied  square 
sheet,  since  all  subsequent  cuts  along  the  parallel,  diagonal,  and 
perpendicular  were  made  relative  to  this  specified  edge.  A  test  was 
performed  with  the  intent  to  verify  the  second  assumption  and  to  align  the 
four  filters  to  their  proper  rotations  relative  to  the  common  reference 
plane. 


The  filter  mount  assembly  was  placed  over  a  diffuse,  unpolarized 
illumination  source.  A  second  HN38  square  sheet  was  mounted  over  a 
protractor  (refer  to  Figure  3.1-1)  with  the  corners  fixed  at  the  four 
ordinates.  The  sheet  was  rotated  until  a  relative  minimum  transmission 
(maximum  absorption)  via  cross-polarization  was  detected  with  a  light 
meter  for  each  of  the  four  filters.  For  each  rotation,  the  angle  on  the 
protractor  was  measured  to  within  0.5  degree.  This  process  was  repeated 
three  times  and  the  measurements  then  averaged  and  examined  for  bias. 
The  standard  error  was  within  the  measurement  accuracy;  no  rotational 
bias  was  detectable  for  any  of  the  four  filters  within  the  measurement 
accuracy  of  this  test. 

Figures  3.1-2  through  5  illustrate  the  alignment  process  for  each  of  the 
four  filters.  These  figures  also  serve  as  a  simple  demonstration  of  the  law 
of  Malus:  transmission  of  light  through  the  two  overlapping  linear  sheet 
polarizers  is  proportional  to  the  square  of  the  cosine  of  the  angle  between 
their  principal  axes  of  transmission. 
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Figure  3.1-4.  Cross-polarized  90°  Figure  3.1-5.  Cross-polarized  135° 

filter  (second  filter  from  bottom).  filter  (top  filter). 


Film  calibration 

One  type  of  film  was  used  for  the  entire  data  collection:  36- 
exposure  Kodak  Gold  200-speed  color  negative  film.  This  film  was  bulk 
purchased  from  a  single  vendor;  and,  the  cartons  were  inspected  to 
ensure  that  the  film  rolls  had  the  same  manufacturing  lot  number. 

A  sensitometric  24-step  wedge  was  exposed  on  the  leading  edge 
of  each  film  roll  with  a  Joyce-Gevaert  type  2L  sensitometer.  A  frame  mask 
overlay  was  created  that  would  simulate  the  frame  placement  of  the 
Nashika  camera  and  yet  include  coverage  of  all  24  exposure  steps.  The 
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intent  of  this  mask  was  to  support  digital  image  transfer  of  the  full  step 
wedge  during  photo  CD  processing. 

3.1.2  In-field  imaging 

Overview 


The  intent  of  in-field  imaging  was  to  collect  polarimetric  data  on  a 
set  of  controlled  test  surfaces  under  natural  illumination  conditions:  in  this 
experiment,  under  polarized  clear  sky  radiance.  The  camera  orientation 
was  fixed  with  respect  to  solar  azimuth,  the  test  surfaces  were  rotated 
through  four  discrete  angles,  and  images  were  collected  at  four  distinct 
solar  elevations. 

Construction  of  the  in-field  test  stand 

A  test  stand  was  constructed  for  the  purpose 
of  fixing  the  imaging  polarimeter  over  a  set  of  test 
surfaces  in  the  field. 

Figure  3.1-6  illustrates  the  full  setup  of  the  in-field  test 
stand  (under  less  than  ideal  illumination  conditions). 

A  4x8  foot  (122x244  cm)  plywood  sheet  Figure  3.1-6.  Setup 

of  the  test  stand. 

served  as  the  base  of  the  test  stand.  Four  10-foot 
(3.05  m)  sections  of  electrical  conduit  formed  the  legs.  A  plastic  flange 
provided  a  rigid  camera  mount  at  the  apex  of  the  test  stand.  The  resulting 
effective  focal  distance  between  the  camera  system  aperture  and  the  test 
stand  base  is  nine  feet  (2.74  m).  This  distance  provided  reasonable 
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mitigation  of  the  parallax  effects  created  by  the  finite  separation  of  the  four 
lenses  (55  mm  between  the  two  most  distant  lens  centers). 

Figure  3.1-7  shows  the  alignment  of  the 
imaging  polarimeter  over  the  in-field  test 
stand. 

Construction  of  the  test  surface  panels 

Two  identical  2x4  foot  (61x122  cm) 
lighting  diffuser  panels  were  used  as  test 
surface  panels.  The  unique  feature  of  these 
panels  is  their  ‘cracked  ice’  pseudo-random  texture,  which  ensures  a  high 
degree  of  spatial  correlation.  To  ensure  that  the  experiment  contained  at 
least  one  surface  that  was  within  the  linear  range  of  the  characteristic 
curves,  one  panel  was  coated  with  chrome  aluminum  spray  paint  and  the 
other  was  coated  with  gray  primer.  The  two  panels  were  mounted  on  a 
4x4  foot  (122x122  cm)  plywood  sheet  with  white  thumbtacks,  which  also 
served  as  in-scene  ground  control  points  during  image  co-registration. 

The  entire  assembly  was  then  mounted  on  the  center  of  the  test  stand 
with  a  single  bolt  through  the  center  of  both  plywood  sheets.  This 
mounting  strategy  allowed  for  360'’  rotation  of  the  test  surface  panels  in 
relation  to  the  fixed  reference  plane  of  the  imaging  polarimeter. 

Other  test  stand  features 

An  18-inch  (45.7  cm)  back-surfaced  acrylic  hemispherical  security 
mirror  was  used  to  reflect  a  polar  representation  of  the  complete  skydome 
at  the  beginning  of  each  collection.  A  20-foot  air-driven  shutter  release 
was  employed  to  minimize  in-scene  obscuration  of  the  sky  radiance  by 


Figure  3.1-7.  Alignment  of 
the  imaging  polarimeter  on 
the  in-field  test  stand. 
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the  photographer.  A  5.5-inch  (14  cm)  vertical  post  was  mounted  on  the 
front  (sun-facing)  edge  of  the  base  in  order  to  cast  a  linear  shadow  that 
could  support  test  stand  alignment  with  the  changing  solar  azimuth. 
Finally,  a  Kodak  10-step  gray  scale,  8-step  color  control  patch,  and 
neutral  density  card  was  included  in  each  image  collection  to  support  in¬ 
scene  secondary  calibration. 

Figure  3.1-8  illustrates  the  two 
base-mounted  test  panels  and  the 
several  test  stand  features. 

Mission  planning  for  the  in-field 
collection 

A  key  requirement  for  optimal  in-field  data  collection  was  to  obtain 
clear  sky  illumination  during  the  collection  period,  nominally  a  four-hour 
window  that  would  either  begin  or  end  at  solar  noon. 

Hourly  visible  and  infrared  weather  satellite  (GOES)  imagery  was 
available  from  a  University  of  Illinois  (Champaign-Urbana)  database 
hosted  on  the  Internet.  Based  on  a  one-week  forecast  of  regional  high 
pressure  fronts  during  the  first  days  of  November  1994,  this  satellite 
imagery  was  downloaded  in  support  of  near-term  prediction  of  an  optimal 
collection  window  for  the  in-field  test.  Weather  images  from  the  early 
morning  of  07  Nov  94  indicated  that  such  a  collection  window  would  exist 
over  Syracuse  NY  for  the  entire  daylight  period.  The  in-field  collection  was 
subsequently  scheduled  to  begin  at  solar  noon,  approximately  1155  EST. 
The  intent  was  to  provide  the  largest  range  of  solar  elevations  beginning 
with  the  maximum  elevation  at  solar  noon  and  the  minimum  at  sunset. 


Figure  3.1-8.  In-field  test  panels  and 
supporting  test  stand  features. 
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Figure  3.1-9.  Visible  satellite  image,  1200  Eastern  Standard  Time. 


Figure  3.1-10.  Infrared  satellite  image,  0800  Eastern  Standard  Time. 
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Figures  3.1-9  and  10  illustrate  the  presence  of  a  high-pressure  front  with 
cloud-free  skies  over  the  Syracuse  NY  area  on  07  November  1994.  (Note 
the  large  tropical  storm  over  the  Atlantic.) 

Execution 


In-field  polarized  photographs  were  taken  on  a  near-cloudless  day, 
07  Nov  94,  near  the  Skytop  parking  lot  on  the  Syracuse  University  NY 
campus.  This  location  afforded  nearly  horizon-to-horizon  visibility  at  most 
azimuths.  Six  sets  of  photographs  were  taken  over  a  six-hour  period, 
beginning  one  hour  before  solar  noon  and  then  at  subsequent  one-hour 
intervals.  The  entire  platform  was  rotated  slightly  for  each  collection  so 
that  the  polarimeter  reference  plane  remained  parallel  with  the  principal 
plane:  a  plane  containing  the  sun  and  a  vertical  line  running  through  the 
base  center,  the  polarimeter,  and  local  zenith.  The  shadow  cast  by  the 
small  vertical  post  mounted  on  the  platform  was  used  to  maintain  this 
alignment.  The  length  of  the  post  shadow  was  also  used  to  provide  an  in¬ 
scene  secondary  measurement  of  solar  elevation. 

Each  of  the  six  datasets  consisted  of  five  collections.  For  the  first 
collection,  the  dome  mirror  was  placed  over  the  center  of  the  test  stand 
base.  The  resulting  collection  provided  a  nearly  complete  image  of  the  full 
skydome  reflecting  off  the  mirror  and  onto  the  four  focal  planes  of  the 
overhead  polarimeter. 
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Figure  3.1-1 1  illustrates  the  first  test  scene, 
containing  a  dome  mirror,  that  was  collected 
for  each  of  the  four  datasets.  Note  the 
substantial  cloud  cover  that  is  apparent  in  the 
skydome  reflection,  indicative  of  an  earlier 
unsuccessful  collection  attempt  that  was  not 
included  in  the  data  analysis  or  results. 

For  the  remaining  four  collections,  the 
test  panels  were  rotated  0°,  45°,  90°,  and 
135°  with  respect  to  the  defined  reference  plane.  The  intent  of  this  rotation 
was  to  vary  the  illumination  geometry  of  the  scene  relative  to  the  test 
target,  an  effect  which  would  normally  be  difficult  (if  not  impossible)  to 
achieve  with  most  natural  scene  contents  in  the  terrestrial  environment 
(e.g.,  a  forest  or  lake). 

Figures  3.1-12  through  15  illustrate  examples  of  the  four  remaining  test 
scenes  that  were  collected  for  each  of  the  four  datasets.  The  test  panels 
were  succesively  rotated  counterclockwise  45°  for  each  collection:  test 
scene  2  is  at  0°  rotation,  scene  3  is  at  135°,  scene  4  is  at  90°,  and  scene 
5  is  at  45°.  These  examples  are  taken  from  the  first  image  in  each  of  the 
quad  series,  i.e.,  the  quad  image  with  the  polarizer  rotated  90°  from  the 
reference  plane. 


Figure  3.1-11.  Test  scene  1. 
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Figure  3,1-14^  Test  scene 

Figure  3.1-15.  Test  st^ne  5 

3.1.3  Aerial  Imaging 

Overview 

The  intent  of  aerial  imaging  was  to  collect  polarimetric  data  under 

clear  skies  of  large  natural  surface  features  containing  resolvable 
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Figure  3.1-12.  Test  scene  2. 

Figure  3.1-13.  Test  scene  3. 

homogeneous  texture,  in  particular,  the  waterwave  surface  of  the  Hinckley 
reservoir  near  Marcy  NY.  Other  features  of  possible  future  research 
interest  would  be  collected  as  targets  of  opportunity.  Because  of  the 
relative  lack  of  control  of  the  collection  geometry,  compared  with  the  in¬ 
field  test,  much  of  the  subsequent  data  analysis  would  qualify  this 
collection  as  a  demonstration. 

Mission  planning  for  the  aerial  collection 

The  opportunity  to  collect  aerial  polarimetric  imagery  presented 
itself  with  a  short  timetable.  The  flight  was  scheduled  for  03  November 
1994  based  on  24-hour  weather  information:  sunny,  lightly  scattered 
clouds,  and  light  winds.  Two  missions  would  be  flown,  one  occuring  close 
to  solar  noon  and  the  other  approximately  two  hours  later. 

Description  of  the  flight  platform 

A  Cessna  Skyhawk  single-engine  airplane  was  rented  from 
Landcare  Aviation,  Inc.,  located  at  the  Marcy  NY  regional  airport.  The 
aircraft  had  two  inspection  ports  in  the  aft  underbelly  of  the  fuselage,  both 
available  for  use  as  camera  viewports.  The  imaging  polarimeter  was 
mounted  over  one  of  the  viewports  and  aligned  with  its  reference  plane 
roughly  parallel  with  the  flightline;  in  lieu  of  a  rigid  mounting  frame,  the 
polarimeter  was  cradled  on  an  exposed  area  of  the  aircraft’s  foam  rubber 
insulation.  The  aircraft’s  on-board  magnetic  compass  and  altimeter 
provided  the  only  indications  of  camera  azimuth  and  elevation. 

Figures  3.1-16  through  18  show  the  aircraft  used  for  the  aerial  data 
collection,  indications  of  the  weather  and  illumination  conditions,  and  the 
location  of  the  polarimeter  in  the  aircraft. 
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Figure  3.1-16.  L-R:  photographer,  flight  platform,  and  pilot. 


Figure  3.1-17.  Interior  view  of 
polarimeter  placement. 


Figure  3.1-18.  Exterior  view  of 
polarimeter  placement. 


Execution 

The  primary  imaging  target  for  both  missions  was  the  water  surface 
of  the  Hinckley  reservoir,  located  approximately  10  miles  NNE  of  the 
Marcy  NY  airport  and  four  miles  ESE  of  Griffiss  Air  Force  Base  near 
Rome  NY.  Three  passes  of  the  reservoir  were  made  at  three  different 
flying  heights:  approximately  1000,  500,  and  150  feet  above  the  water 
surface.  During  the  return  leg  of  both  flights,  targets  of  opportunity  were 
also  collected:  farm  fields,  tree  stands,  rural  buildings,  et  cetera. 
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Figure  3.1-19.  Subset  of  the  USGS  quadrangle  containing  the  Hinckley  reservoir. 


For  the  first  mission,  10  aerial  images  were  collected  during  the 
interval  1312  -1327:  two  images  were  collected  for  each  of  the  three 
overflights  of  the  reservoir  and  four  images  were  collected  of  opportunistic 
targets.  For  the  second  mission,  10  aerial  images  were  collected  during 
the  interval  1530  -  1548:  again,  two  images  were  collected  for  each  of  the 
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three  overflights  of  the  reservoir  and  another  four  images  were  collected 
of  opportunistic  targets. 


Figure  3.1-20.  Aerial  view  of  the  Hinckley 
reservoir. 
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3.2  Image  Data  Processing 


Film  development  and  image  digitization 

The  film  was  turned  over  to  a  commercial  vendor  for  standard 
KODAK  film  development  and  direct  Photo  CD™  processing.  Special 
processing  instructions  included  requests  to  turn  off  the  automatic  color 
balance  control,  to  scan  in  the  step  wedge,  and  to  center  the  two  pairs  of 
images  within  two  standard  35-mm  frames  -  since  the  dimensions  of  each 
photo  pair  are  approximately  22  mm  x  37  mm.  The  standard  Photo  CD™ 
process  converts  a  normal  35  mm  frame  (24  mm  x  36  mm)  into  a  2048  x 
3076  pixel  image  at  a  scanning  resolution  of  approximately  11.7  microns  x 
11.7  microns.  As  a  result,  one  millimeter  of  exposed  film  from  each  frame 
was  not  scanned,  but  two  millimeters  of  unexposed  film  were  scanned. 

Figures  3.2-1  and  2  illustrate  a  complete  photo  ‘quad’  for  one  data 
collection.  Note  that  each  cross-polarized  pair  exists  as  two  subimages 
that  occupy  the  same  35-mm  image  frame;  the  two  35-mm  frames  that 
form  each  quad  are  sequential  on  the  roll  of  film.  Also  note  the  two  mm  of 
unexposed  film  on  the  upper  edge  of  the  frame  and  the  unexposed 
vertical  band  that  separates  the  two  subimages;  the  one  mm  of  exposed 
but  unscanned  film  is,  for  obvious  reasons,  missing  from  the  left  and  right 
edges  of  the  digitized  image  frames. 


polarized  image  pair.  polarized  image  pair. 
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Creation  of  the  digital  image  quads 


Each  Photo  CD™  compact  disk  was  mounted  on  a  Compudyne 
486  DX2  computer  system  operating  in  a  Windows™  environment.  The 
commercial  software  program  ADAPTEC  Magic  Lantern™  was  used  to 
ingest  each  image  pair,  crop  the  pair  to  create  two  separate  digital 
images,  and  then  convert  each  image  into  a  Tagged  Image  Format  (TIF) 
file  for  transfer  to  a  1-gigabyte  hard  disk  drive.  Each  dataset,  comprising 
four  uncompressed  3-color  TIF  files,  occupied  almost  38  megabytes  of 
hard  disk  space;  the  80  TIF  files  generated  by  the  in-field  collection 
required  almost  756  megabytes. 

The  software  program  NOVASTOR  Novatar™  was  used  to  copy  all 
TIF  files  from  the  Compudyne  hard  disk  drive  to  8-mm  data  tape.  The  files 
were  put  into  a  TAR  format  acceptable  to  a  SUN  UNIX™  environment. 

Co-registration  of  the  digital  image  quads 

The  software  program  SRI  Environment  for  Visualizing  Images 
(ENVI™)  was  used  for  the  remainder  of  the  data  processing.  It  was 
hosted  on  a  SUN  Sparc™  2  workstation  operating  in  a  UNIX™ 
environment.  All  computer  systems  and  software  utilized  in  this  study  are 
part  of  the  FIRS^  cluster  located  in  Bray  Hall  on  the  SUNY  ESF  campus. 


^  Facility  for  Image  Processing  &  Remote  Sensing  at  SUNY  ESF 
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Figure  3.2-3 
illustrates  the 
SUN  workstation 
configuration  that 
was  utilized  for 
the  majority  of  the 
digital  image 
processing. 


Figure  3.2-3.  SUN  workstation  and  peripherals  at  the  FIRS. 


The  four  images  from  each  collection  were  co-registered,  using  in¬ 
scene  markers  (thumbtacks)  as  common  ground  control  points  and  then 
selecting  a  first-order,  four-point  (rotation,  scaling,  &  translation)  solution 
within  ENVI.  Due  to  the  close  imaging  geometries  of  the  four  images,  co¬ 
registration  was  predominantly  a  simple  translation  of  image  coordinates: 
standard  errors  of  less  than  0.30  pixel  were  readily  achievable. 


Creation  of  the  Stokes  images  and  Stokes  derivative  images 

Using  the  ENVI  Band  Math  option,  the  co-registered  images  were 
mathematically  combined  to  form  the  Stokes  parameter  images: 


I  =  [lmg(0°)  +  Img  (45°)  +  lmg(90°)  +  lmg(135°)]  /  2 

Q  =  [lmg(0°)  -  lmg(90°)] 

U=  [lmg(45°)-lmg(135°)] 


Equation  3-1 


Equation  3-2 


Equation  3-3 
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The  Stokes  parameter  images  were  then  mathematically  combined 
to  form  the  Stokes  derivative  images,  P,  the  percentage  polarization,  and 
T,  the  orientation  angle  of  the  polarization  ellipse: 


p  = 

J 

*100[%]  ,and 

Equation  3-4 

V  "  J 

/  \ 

T  = 

u  it/ 
—  *  tan  *  — 

*  —  [deg] 

Equation  3-5 

\2  Qy 

n 

Also,  since  the  sum  of  each  orthogonally  polarized  image  pair  represents 
an  estimate  of  the  unpolarized  intensity  image  (I),  the  histogram  statistics 
of  the  difference  of  the  two  sums  were  calculated  in  order  to  provide  a 
system  calibration  check.  If  the  system  was  perfectly  calibrated,  the 
difference  image,  D,  should  be  everywhere  zero: 


D  = 


r  - 1”  =  [lmg(0°)  +  lmg(90°)]  -  [lmg(45°)  +  lmg(135°)] 


Equation  3-6 


Other  derivative  images  were  also  calculated: 


IN  =  the  difference  between  the  maximum  image  digital  count  and  I 
(effectively,  a  digital  negative), 

iM-  .  I  Equation  3-7 


PN  =  the  percent  depolarization, 
PN=  100 -P[%] 


Equation  3-8 
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TC  =  the  90°  complement  of  T, 


Equation  3-9 


Also,  the  45°  and  135°  complements  of  T  were  calculated  for  subsets  of 
the  in-field  collection  images: 

TB  =  the  45°  complement  of  T, 


TB  =  —  * 
n 


tan 


-I 


tein 


(r+45‘’)/  —  ,,,, 
f  K180JJJJ 


[deg] 


Equation  3-10 


and 

TD  =  the  135°  complement  of  T, 


* 

1 

00 

ii 

Q 

1- 

r 

tan"^ 

tanf(r  +  135'')/f^ll 

[deg] 

Equation  3-11 

TC 

1 

y 

Based  on  considerations  elaborated  in  Chapters  1  and  2,  only  the 
red  color  band  of  the  three  spectral  bands  was  processed  and  retained  for 
image  data  analysis.  One  additional  consideration  that  came  to  light 
during  image  processing  was  the  data  limit  imposed  by  the  64-megabyte 
capacity  of  the  workstation’s  random  access  memory  (RAM).  The  first  six 
equations  in  the  preceeding  section  were  used  to  process  each  full  co¬ 
registered  image  quad,  creating  an  effective  six-image  real  array  with 
nominal  dimensions: 
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6  images  x  (1400  x  1850)  pixels/image  x  4  bytes/pixel  =  62,160,000  bytes. 


However,  this  array  represents  a  single  color  band.  While  ENVI  provides 
for  several  options  which  store  image  bands  on  disk  versus  RAM,  it 
became  apparent  that  the  logistics  of  creating  one-band  image  files  and 
then  aggregating  multiple  one-band  image  files  into  single  multiband 
image  files  would  be  impractical  for  large  numbers  of  three-color-band 
quad  images  -  effectively,  images  containing  12  co-registered  bands. 
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3.3  Image  Data  Simulation 


A  synthetic  imaging  model  was  created  to  simulate  Stokes  imagery 
of  the  reflected,  polarized  skydome  under  clear  sky  conditions  to  compare 
with  the  first  scene  of  each  of  the  four  datasets  of  the  in-field  test  (refer  to 
Figure  3.1-1 1).  This  section  describes  the  composition  of  the  image  data 
simulation. 


In  addition,  the  FORTRAN  source  code  that  was  developed  in  this  study 
for  the  simulation  is  listed  in  Appendix  A.  The  code  is  written  in  Microsoft 
FORTRAN  for  Windows™;  the  executable  files  run  under  MS  Windows™. 

Calculation  of  clear  sky  radiance 

A  typical  clear  skydome  radiance  distribution  is  the  CIE  model  [CIE, 

1973]: 


where 

^  =  0.91 +  10e"''‘+ 0.45  cos  V 

_ Y _ ^”0.32/cos6 

and 

C  =  0.274(0.91  +  lOe"^®" )  +  0.45  cos^  0  „ 


Equation  3-12 


Equation  3-13 


Equation  3-14 


Equation  3-15 
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such  that 


e  =  the  zenith  angle  to  a  sky  element, 
cp  =  the  azimuth  angle  to  a  sky  element, 

II  =  the  angle  from  the  sun  center  to  a  sky  element, 

H  =  cos“’[cos0o  COS0  +sin0o  sin0  cos((p  -(t)o)] 


Equation  3-16 


00  =  the  zenith  angle  to  the  sun  center,  and 
(po  =  the  azimuth  angle  to  the  sun  center. 

Although  the  CIE  model  provides  broadband  spectral  radiance  values, 
results  from  this  model  must  be  used  with  caution  when  applied  to 
narrowband  spectral  radiance,  particularly  in  the  vicinity  of  the  solar 
spectral  radiance  peak  at  460  nanometers. 

Calculation  of  clear  sky  polarization 

The  direct  radiance  from  the  sun  is  fully  unpolarized.  However,  the 
fraction  of  unpolarized  solar  radiance  that  becomes  partially  polarized  by 
atmospheric  Rayleigh  scattering  can  be  described  by  the  function; 


0.94 

2  > 

1-COS  la 

2 

Equation  3-17 

U  +  cos  ixj 

This  function  accounts  for  multiplicative  depolarization  created  by  the 
anisotropy  of  atmospheric  scatterers  [Coulson,1975]:  in  this  model, 
fractional  polarization  varies  from  0  (unpolarized)  in  the  vicinity  of  the  sun 
to  .94  (94%  polarized)  at  the  anti-solar  point. 
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Calculation  of  the  azimuth  angle  of  the  polarization  ellipse 


The  direction  of  polarization  is  always  normal  to  a  plane  that 
contains  both  the  incident  radiance  vector  from  the  sky  element  of  interest 
and  the  incident  solar  radiance  vector.  These  two  vectors  are  first 
analyzed  into  their  unit  vector  components: 

Q  =  the  sky  incident  radiance  vector 


a 

sinG  coscp 

Q  = 

e. 

= 

sinG  sincp 

Equation  3-18 

Q.. 

cosG 

S  =  the  solar  incident  radiance  vector 


's; 

sinGo  cos(Pq 

S  = 

sinG  0  sin  (po 

Equation  3-19 

A. 

cos  Go 

The  polarization  vector  P  is  the  cross  product  of  S  and  Q: 


P  = 

Py 

=  SxQ  = 

S.Q.-S^Q, 

S.Qy-S^Q, 

Equation  3-20 


This  cross-product  specifies  the  orientation  of  the  polarization  ellipse  for 
the  pure  Rayleigh  scattering  model. 

For  the  in-field  experiment,  the  reference  axis  is  the  z-axis  and  the 
reference  plane  is  the  y-z  plane  (the  sun  is  always  in  the  +y  direction). 
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Therefore,  the  calculation  of  the  azimuth  angle  of  the  polarization  ellipse 
for  each  sky  element  is: 


®poi  =  tan 


-i 


\PyJ 


Equation  3-21 


If  a  sky  element  is  within  the  reference  plane  (which  also  coincides  with 
the  principal  plane  for  the  in-field  test),  the  azimuth  angle  is  0°. 

Calculation  of  the  incident  radiance  Stokes  vector 


The  incident  radiance  Stokes  vector  Sgky  for  each  sky  element  is 
readily  calculated  using  Equation  2-28: 


9)  = 


■  I ' 

L 

Q 

IPcos2©p„, 

u 

ZPsin20^„, 

_«  0 

«0 

Equation  3-22 


Calculation  of  the  reflected  Stokes  vector 

The  Stokes  vector  image  resulting  from  this  imaging  system  can  be 
fully  described  by  Mueller  matrix  calculus: 


Sou,  =  [Polarizer{Q)\Rot{a)\Mirror{^)\Rot{-a)\S^^^ 


Equation  3-23 


where 

Rot(-a)=  rotation  of  the  incident  Stokes  vector  Sgky  from  the  polarizer 
reference  frame  to  the  mirror  reference  frame  (the  azimuth  anglea  of 
each  facet). 
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0 


0 


Rot  (-a) 


1  0 
0  cos(-a)  sin(-a)  0 
0  -sin(-a)  cos(-a)  0 
0  0  0  1 


Equation  3-24 


Mirrori^)  =  the  polarizing  effect  of  reflection  from  the  mirror, 


Mr(|5)  =  | 


0 

0 


0 

0 


0 

0 

-2r^n 


0 

0 

0 

-^Vl 


Equation  3-25 


Rot(a)  =  rotation  of  the  mirror-reflected  Stokes  vector  back  to  to  the 
polarizer  reference  frame, 


Rot(+a)  = 


10  0  0 
0  cos(+a)  sin(+a)  0 
0  -sin(+a)  cos(+a)  0 


0  0  0  1 


Equation  3-26 


and 

Polarizer (Q)  =  the  polarizing  effect  of  each  linear  polarizer. 
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Equation  3-27 
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A  detailed  description  of  the  underlying  mathematics  can  be  found  in 
Clarke  and  Grainger  [1971].  Although  the  parameter  V  (which  quantifies 
the  fraction  of  right-handed  circular  polarization)  is  negligible  in  most 
terrestrial  applications,  it  is  recognized  that  the  aluminum  coating  on  the 
mirror  might  induce  a  small  amount  of  circular  polarization.  However,  the 
fact  that  the  mirror  is  back-surfaced  on  an  acrylic  matrix  suggests  that  any 
otherwise  measurable  amount  of  circular  polarization  will  be  masked 
through  depolarization;  the  results  bear  out  this  hypothesis. 
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3.4  Image  Data  Analysis 


This  section  describes  the  analytic  techniques  that  were  applied  to 
the  Stokes  imagery  obtained  from  both  the  in-field  and  aerial  image  data 
collections.  The  majority  of  the  image  data  analysis  was  accomplished 
with  the  use  of  ENVI™  on-line  analytic  software  capabilities. 

3.4.1  In-field  image  data  anaiysis 

Full  image  scene  analysis 

Means,  variances,  and  histograms  of  pixel  intensity  were  calculated 
for  full  scenes  from  the  calibrated  image  datasets.  These  statistics  provide 
coarse  benchmarks  for  comparison  with  the  sub-scene  (feature)  statistics 
that  follow. 

Dome  mirror  scene  analysis 

512  X  512  pixel  subsets  of  the  dome  mirror  subimages  in  the  Scene 
1  images  were  cropped  for  qualitative  comparison  with  the  synthetic 
imagery  from  the  image  data  simulation. 

In-scene  step  wedge  and  neutral  density  card  analysis 

Means  and  variances  of  pixel  intensity  were  calculated  for  10  x  10 
pixel  subsets  from  each  of  the  10  steps  on  the  in-scene  step  wedge  that  is 
visible  in  each  collection.  Means  and  variances  of  pixel  intensity  were 
calculated  for  66  x  66  pixel  subsets  from  the  neutral  density  card  visible  in 
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each  collection.  Figure  3.4-1  illustrates  the  ENVI™  capability  for  selecting 
small  regions  for  statistical  analysis. 


Figure  3.4-1.  In-scene  step 
wedge  and  neutral  density  card 
measurement. 


Figure  3.4-2.  In-scene  solar 
elevation  measurement. 


In-scene  solar  elevation  measurement 


The  length  of  the  shadow  formed  by  a  small  (5.5-inch)  vertical  post 
mounted  on  the  platform  was  measured  using  the  ENVI™  cursor  location 
capability.  The  pixel  length  was  calibrated  with  the  known  length  of  the 
panels  in  order  to  accurately  estimate  the  solar  elevation.  Figure  3.4-2 
illustrates  the  visible  shadow  in  a  negative  intensity  image. 

The  solar  elevation  is  then  directly  calculated: 


Solar  Elevation  =  tan 


^Shadow  Length[m^ 
5.5  [in]  y 


Equation  3-28 


Panel  scene  analysis 

Means,  variances,  and  histograms  of  pixel  intensity  were  calculated 
for  large  representative  areas  of  the  textured  panels.  The  statistics  were 
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collected  for  both  fully  illuminated  and  mirror-shadowed  areas  in  the 
Scene  1  images. 

Power  Spectral  Density  (PSD^  analysis 

Representative  128  x  128  pixel  subsets  of  the  two  panels  were 
cropped  from  Scenes  2  through  6  for  PSD  analysis.  The  means, 
variances,  and  histograms  were  first  calculated  from  the  image  subsets 
prior  to  fast  Fourier  transformation  (FFT)  and  generation  of  the  PSDs. 

In  order  to  establish  a  basis  for  comparison  between  scenes,  normalized 
power  spectral  densities  (PSDs)  were  calculated  for  each  subimage: 


PSD„„„„{K^ky) 


[FW, /)J 

Equation  3-29 

2  X  l\j\J /O 

S 

The  mean  intensity  is  first  subtracted  from  the  image.  The  motivation  for 
subtracting  the  mean  is  to  make  the  central  ordinate  of  the  FFT  zero. 
Since  the  central  ordinate  of  the  FFT  is  the  mean,  does  not  contribute  to 
the  evaluation  of  image  variance,  and  generally  has  the  greatest 
magnitude  of  all  the  spectral  components,  this  step  facilitates  scaling  of 
the  transformed  image: 


FFT, ^(0,0)  =  7  =  0 


Equation  3-30 


The  unbiased  image  is  then  2D  fast  Fourier  transformed  through  the 
ENVI™  filter  option.  The  PSD  is  created  by  squaring  the  transform: 


PSD{k„k,)  =  \FFT{K,k,i 


Equation  3-31 
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The  third  step  divides  the  PSD  by  the  intensity  variance  of  the  image.  The 
motivation  is  to  scale  the  PSD  so  that  each  spectral  component  now 
represents  its  fractional  contribution  to  the  total  variance,  since 


\\PSD{k^,k^)dk,dk^=s^ 


Scaling  the  PSD  by  100%  converts  the  fractional  contribution  of  each 
spectral  component  to  a  percent  contribution. 

PSD  stability  analysis 

A  primary  objective  of  this  study  is  to  compare  the  relative  stability 
of  polarimetric  power  spectral  estimation  over  the  given  range  of  imaging 
geometries.  In  the  absence  of  more  robust  stability  analysis  techniques, 
existing  ENVI™  capabilities  provide  for  an  averaged  periodogram 
approach,  as  described  by  Kay  [1988].  The  underlying  assumptions  in 
applying  this  approach  are  1)  the  PSD  from  each  realization  represents  an 
estimate  of  a  similar  process,  in  this  case,  reflected  intensity  from  the 
same  pseudo-random  (and  therefore  deterministic)  surface,  i.e.,  the 
imaged  scenes  of  the  ‘cracked  ice’  diffuser  panels;  and  2)  significant 
deviations  from  the  averaged  PSD  are  the  result  of  differences  in  imaging 
geometry. 

In  order  to  establish  a  basis  of  comparison  between  image  spectra,  the 
images  were  first  rotated  to  a  common  azimuth  with  respect  to  panel 
orientation,  taking  into  account  the  bias  error  caused  by  rotation  of  the 
polarimeter.  128  x  128  pixel  subimages  of  the  panels  were  cropped  and 
the  statistics  (mean  and  variance)  calculated.  As  before,  the  mean  was 
subtracted  from  each  image  prior  to  calculation  of  the  PSD  and  the  PSDs 


Equation  3-32 
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normalized.  An  average  PSD  was  calculated  for  each  combination  of 
imaging  geometry  factors,  azimuth  angle  and  solar  elevation  angle: 

PSD  (eh  A)  —  ~7^^PS^(ehA,azB) 

^  B=\ 


PSD(azB)  —  PSD(gh  A.az  B'S 
^  A=\ 

The  Averaged  Difference  PSD  (ADP)  was  then  calculated  for  each 
combination  of  factors: 


Equation  3-34 


Equation  3-33 


ADR 


{eh  A) 


4^1  ( 


{eh  A,az  B) 


5=1 


-  PSD(ehA)^ 


Equation  3-35 


'iDP^.„=\t\PSD, 

^  A=\ 


{elvA.azB) 


—  PSD{azB)\ 


Equation  3-36 


This  spectrum  represents  the  average  difference  in  the  variance  at  each 
spectral  component  between  the  averaged  PSD  and  each  of  the  PSDs 
that  contributed  to  the  average;  it  represents  a  relatively  straightforward 

estimate  of  the  stability  of  PSD .  The  figure-of-merit  for  PSD  stability  is  the 
integrated  difference  variance  (IDV): 


IDV  -  j^ADP(k^,ky)dk^dky 


In  lieu  of  a  packaged  ENVI™  capability  to  execute  a  double  integration 
over  a  scalar  ADP  array,  an  alternative  approach  was  to  run  the  ENVI™ 
statistics  tool,  create  a  histogram  of  the  ADP,  transfer  the  histogram  data 


Equation  3-37 
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to  a  file,  and  then  run  a  simple  FORTRAN  program  to  integrate  over  the 
histogram: 


IDV  =  '^(Pixel  Count)(AVar^) 

w=0 


Equation  3-38 


The  deficit  of  using  this  approach  is  that  the  histogram  bins  AFar  into 
discrete  values,  creating  an  unknown  approximation  error  for  values  in  the 
lowest  bin.  Since  the  expectation  is  that  a  majority  of  values  are  near¬ 
zero,  all  the  variances  below  the  value  of  the  next  lowest  bin  are  truncated 
to  zero  by  the  histogram  and  there  is  no  way  to  evaluate  them.  However, 
there  is  a  related  side  benefit  to  this  approach:  on  the  assumption  that  it  is 
the  larger  values  of  AVar  that  indicate  significant  trends,  a  raised  (non¬ 
zero)  threshold  can  be  used  to  evaluate  a  fractional  IDV: 


IDV^racfionai  =  S  Count)(AVarJ 

w- threshold 


Equation  3-39 


Several  threshold  values  are  used  in  order  to  demonstrate  the  validity  of 
this  histogram  approach  in  the  calculation  of  IDV. 


3.4.2  Aerial  image  data  analysis 

While  leaving  the  details  to  Chapter  4,  this  section  must  address 
the  fact  that  only  one  of  the  12  over-water  collections  was  fully  processed 
and  available  for  polarimetric  image  analysis.  The  results  from  this 
truncated  analysis  of  the  one  aerial  image  are  provided  as  a  qualified 
demonstration  of  feasibility  rather  than  as  an  experimental  outcome. 
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Homeogeneous  sub-scene  analysis 


Means,  variances,  and  histograms  of  pixel  intensity  were  calculated 
for  one  large  (512  x  512  pixel)  representative  area  (wind-driven  water- 
wave  surface)  that  contained  homogeneous  spatial  content. 

Power  Spectral  Density  fPSD)  analysis 

Four  256  x  256  pixel  subsets  were  cropped  from  the  larger  sub¬ 
scene  for  PSD  analysis.  The  means,  variances,  and  histograms  were  first 
calculated  from  the  image  subsets  prior  to  Fast  Fourier  Transformation 
(FFT)  and  generation  of  the  PSDs.  The  normalized  PSDs  were  then 
generated  using  the  same  method  described  for  the  in-field  analysis. 

PSD  stability  analysis 

In  the  absence  of  experimental  data  relating  the  imaging  geometry  factors 
of  elevation  and  flight  azimuth,  an  average  PSD  was  calculated  for  the 
four  spatially  contiguous  sub-scene  PSDs: 

PSD(suh-scene)  =  —  PSD(suhsceneC'i 
^  C=1 


This  analysis  should  demonstrate  the  assumption  of  spatial  invariance  for 
the  quasi-stationary  pseudo-Gaussian  process  of  a  localized  wind-driven 
waterwave  surface  [Kinsman,  1965].  This  analysis  should  also  provide  yet 
a  third  averaged  periodogram  method  for  comparing  the  relative  PSD 
stabilities  of  the  polarimetric  image  components. 


Equation  3-40 
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4.0  Results  and  discussion 


The  full  description  of  results  and  discussion  comprises  five 
sections: 

General  results  -  This  section  describes  the  analysis  of  the  general 
execution  of  the  collection  and  processing  phases. 

In-field  calibration  results  -  This  section  describes  the  calibration 
analysis  of  the  in-field  imaging  experiment. 

Image  simulation  results  -  This  section  describes  the  qualitative 
comparison  of  the  imaging  simulation  with  real  imagery. 

Aerial  calibration  results  -  This  section  describes  the  calibration 
analysis  of  the  aerial  imaging  demonstration. 

Spatial  analysis  results  -  This  section  describes  the  spatial  analysis 
of  sub-image  features  in  both  the  in-field  and  aerial  collections. 


4.1  General  results 

Quad  image  co-registration 

The  four  members  of  each  image  quad  were  co-registered  using 
the  registration  tools  available  within  ENVI™.  For  the  in-field  images,  a 
black  dot  in  the  center  of  each  of  eight  in-scene  markers  (thumbtacks) 
served  as  a  common  ground  control  point  in  each  of  the  quad  members 
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For  the  aerial  images,  spatially  distributed  high-contrast  corners  and 
edges  were  used  to  the  maximum  possible  extent. 

Pre-registration  evaluation  of  the  marker  pixel  coordinates  for  all  20 
in-field  image  quads  indicated  that  a  simple  translation  of  image  pixel 
coordinates  was  the  only  requirement  for  co-registration.  This  requirement 
was  consistent  with  the  expectations  for  the  specified  imaging  geometry: 

1)  the  exterior  orientation  aligns 
the  four  optic  axes  parallel  with  nadir 
and  aligns  the  web  axis  parallel  with 
the  defined  reference  plane,  i.e.,  the 
orientation  angles  omega(co),  phi  (cp), 
and  kappa  (k)  are  zero  within  small 
tolerances:  and  rotation  errors  (that 
were  detected  during  analysis)  that 
make  kappa  non-zero  would  cause  a 
uniform  change  in  translation  to  all 
four  quad  members  during  quad  co¬ 
registration, 

2)  the  primary  imaging  target  (the  textured  flat  panels)  were 
centered  on  the  nadir  point  and  aligned  perpendicular  with  the  optic  axes, 

and 


Figure  4.1-1.  Exterior  orientation  coor¬ 
dinate  system  [from  Moffitt  &  Mikhail, 
19801. 


3)  the  four  image  quad  members  were  exposed  on  the  same  spool 
of  film;  the  single  film  spool  is  processed  in-line  on  the  same  digitizer; 
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centering  errors  would  only  cause  a  change  in  x-axis  translation  between 
frames  (image  pairs)  during  quad  registration. 

The  overall  root-mean-square  (RMS)  error  for  the  60  co¬ 
registrations  required  for  the  20  in-field  quad  images  was  0.29  pixel  with  a 
standard  deviation  of  0.34  pixel;  the  overall  RMS  error  for  the  six  co¬ 
registrations  required  for  the  two  aerial  quad  images  (of  the  nine  that  were 
successfully  transferred  to  photo  CD)  was  0.45  pixel  with  a  standard 
deviation  of  0.37  pixel.  This  result  was,  for  the  most  part,  due  to  the  lack 
of  sharp,  high-contrast  common  ground  points.  The  extent  to  which  low 
resolution  optics,  atmospheric  attenuation,  and  aircraft  motion  were 
underlying  causes  is  currently  unknown  but  worthy  of  future  analysis. 

The  main  observation  made  during  the  co-registration  process  was  that 
operator  error  must  be  considered  within  the  registration  error  budget.  The 
visual  selection  of  the  centroids  of  the  in-field  markers  and  the  common 
intersections  of  the  ground  features  in  the  aerial  images  was  subjective. 
Considering  that  co-registration  of  polarimetric  imagery  is  predominantly  a 
translation  of  coordinates,  the  operator  could  be  removed  as  a  source  of 
error  by  the  introduction  of  a  semi-automated  local  cross-correlation  tool 
to  detect  the  centroid  coordinates  of  a  common  in-scene  feature  between 
two  images. 

This  topic  will  be  further  developed  in  Chapter  5. 

Primary  characteristic  curve  calibration 

The  intended  primary  method  for  providing  radiometric  calibration 
of  the  imagery  was  to  use  the  sensitometric  24-step  density  wedge  that 
was  exposed  on  the  leading  edge  of  each  film  roll.  The  step  wedge  would 
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be  transferred  to  photo  CD  and  the  digital  number  (0  to  255)  mean  and 
variance  corresponding  to  the  film  exposure  at  each  wedge  step  could 
then  be  directly  calculated  from  the  digital  image  using  the  ENVI  region 
statistics  option. 

Two  errors  in  processing  the  sensitometric  step  wedge  occurred. 
The  automatic  color/contrast  balance  controls  were  not  disengaged  (as 
requested)  during  photo  CD  processing,  resulting  in  incomplete  transfer  of 
the  full  step  wedge  (due  to  automatic  rejection  of  low-contrast  images) 
and  contrast  distortion  (due  to  automatic  contrast  enhancement)  of  the 
few  step  wedge  images  that  were  transferred.  These  errors  rendered  the 
digitized  step  wedges  unuseable  for  direct  digital  radiometric  calibration. 

However,  the  secondary  in-scene  step  wedges  were  available  on 
each  image  for  the  in-field  imaging  experiment  and  were  subsequently 
used  for  backup  characteristic  curve  calibration.  This  backup  calibration 
source  was  crucial  in  the  presence  of  automatic  color/contrast  balance 
control  since  each  image  frame  (containing  a  quad  pair)  was  contrast- 
adjusted  based  on  its  scene  content.  This  loss  of  calibration  was  most 
significant  for  the  aerial  collections  since  the  inclusion  of  secondary  in¬ 
scene  step  wedges  was  not  feasible  for  the  aerial  imagery. 

Base  upon  this  outcome,  there  are  three  options  to  be  considered 
for  future  polarimetric  image  data  processing: 

1)  include  a  calibrated  in-scene  (or  rather,  in-frame)  step  wedge 
within  each  image  frame  on  the  film  so  that  sensitometric  calibration  can 
be  achieved  independent  the  effects  of  any  post-development  processing 
that  may  occur. 


82 


2)  include  a  single  calibrated  step  wedge  on  each  roll  of  film,  as 
was  currently  done,  but  maintain  full  control  of  post-development  film 
processing,  e.g.,  operate  an  in-house  film  scanner/digitizer,  so  that  all 
post-development  effects  are  known  and  controllable, 

or 


3)  eliminate  the  use  of  film  as  an  intermediate  image  data  transfer 
medium  and  instead  use  a  calibrated  focal  plane  detector  array  to  directly 
generate  a  calibrated  digital  image. 

This  topic  will  be  further  developed  in  Chapter  5. 


83 


4.2  In-field  calibration  results 


General  collection  results 


Twenty-eight  photo  quads  (56  35-mm  frames)  were  collected 
between  1048  and  1554  EST  on  07  Nov  1994.  A  summary  of  the 
collection  results  is  listed  in  Table  4-1. 

The  -1st  solar  elevation  collection  attempt  was  intended  to 
duplicate  the  +ist  solar  elevation  attempt  since  the  two  sets  would  be 
taken  at  the  same  approximate  solar  elevation,  the  -1st  set  taken  one 
hour  before  solar  noon  (maximum  solar  elevation)  and  the  +1st  set  taken 
one  hour  after  solar  noon.  However,  the  -1st  collection  experienced 
several  hardware  failures  with  the  camera  mechanism  and  was  rejected: 
the  resulting  images  from  this  collection  window  were  either  double- 
exposed  and/or  contained  the  photographer  as  an  undesirable  in-scene 
feature  during  the  attempt  to  service  the  mechanism. 

The  0th  solar  elevation  collection  attempt  (circa  solar  noon)  had 
one  additional  mechanical  failure  while  imaging  the  first  scene  but  all 
subsequent  imaging  was  successful.  However,  the  full  0th  collection  was 
also  rejected  for  further  processing  based  on  three  considerations. 

First,  the  collection  was  incomplete  with  the  absence  of  the  first  scene; 
second,  the  remaining  collections  were  all  captured  on  one  roll  of  film, 
eliminating  the  need  for  later  inter-roll  calibration  and  analysis;  and  third, 
the  difference  in  solar  elevation  between  the  0th  and  1st  datasets  was 
small  relative  to  the  elevation  differences  between  the  subsequent  sets. 
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Time 

EST 

Solar 

Elev 

Quad 

Nbr 

Table  4-1 

Collection  Results  Summary 
Frame  Scene 

Nbrs  Nbr 

Remarks 

1048 

-1st 

1 

01/02 

1 

shutter  failure 

1050 

-1st 

2 

03/04 

1 

shutter  failure 

1052 

-1st 

3 

05/06 

2 

advance  failure 

1150 

0th 

4 

07/08 

1 

shutter  failure  & 

advance  failure 

1153 

0th 

5 

09/10 

2 

good  take 

1155 

0th 

6 

11/12 

3 

good  take 

1157 

0th 

7 

13/14 

4 

good  take 

1159 

0th 

8 

15/16 

5 

good  take 

1248 

1st 

9 

17/18 

1 

good  take 

1250 

1st 

10 

19/20 

2 

good  take 

1252 

1st 

11 

21/22 

3 

good  take 

1253 

1st 

12 

23/24 

4 

good  take 

1254 

1st 

13 

25/26 

5 

good  take 

1348 

2nd 

14 

27/28 

1 

good  take 

1349 

2nd 

15 

29/30 

2 

good  take 

1351 

2nd 

16 

31/32 

3 

good  take 

1353 

2nd 

17 

33/34 

4 

good  take 

1354 

2nd 

18 

35/36 

5 

good  take 

1448 

3rd 

19 

37/38 

1 

good  take 

1449 

3rd 

21 

39/40 

2 

good  take 

1451 

3rd 

22 

41/42 

3 

good  take 

1453 

3rd 

23 

43/44 

4 

good  take 

1454 

3rd 

24 

45/46 

5 

good  take 

1548 

4th 

19 

47/48 

1 

good  take 

1549 

4th 

21 

49/50 

2 

good  take 

1550 

4th 

22 

51/52 

3 

good  take 

1552 

4th 

23 

53/54 

4 

good  take 

1554 

4th 

24 

55/56 

5 

good  take 

1  _ i 

In  sum,  because  of  the  intensive  processing  requirements,  these 
considerations  favored  the  selection  of  the  remaining  four  collection  sets 
(20  quads  =  40  frames)  for  follow-on  processing  and  analysis. 
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All  the  exposed  film  from  the  in-field  experiment  was  turned  over  to 
a  comercial  vendor  for  standard  KODAK  film  development  and  direct 
Photo  CD™  processing.  Special  photo  CD  processing  instructions 
included  requests  to  turn  off  the  color/contrast  balance  controls,  to  scan  in 
the  step  wedge,  and  to  center  the  two  pairs  of  images  within  two  standard 
35-mm  frames.  Subsequent  analysis  will  show  that  the  balance  controls 
were  not  disabled  during  processing.  One  of  the  direct  consequences  of 
retaining  this  processing  algorithm  was  rejection  of  most  of  the  step 
wedge  images  and  contrast  distortion  of  the  remainder,  i.e.,  secondary 
characteristic  curve  calibration  was  possible  only  for  the  in-field  imagery. 

Upon  receipt  of  the  processed  CDs,  all  subsequent  digital  data 
processing  and  analysis  proceeded  as  described  in  Chapters  3.2  and  3.3. 

Secondary  characteristic  curve  calibration 

The  backup  method  for  providing  radiometric  calibration  of  the  in¬ 
field  imagery  was  to  use  the  in-scene  1 0-step  density  wedge  that  was 
visible  within  each  image  collection.  Means  and  variances  of  pixel 
intensity  were  calculated  for  10  x  10  pixel  subsets  from  each  of  the  10 
steps,  using  the  ENVI  region  statistics  capability. 

Figure  4.2-2  plots  all  20  characteristic  curves  for  the  uncalibrated  intensity 
(I)  images.  This  plot  demonstrates  the  effect  of  contrast  optimization  that 
was  applied  by  the  photo  CD  processing  algorithm. 
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Scenes  3  and  5  projected  less  total  reflected  radiance  to  the  sensor 
because  the  mean  azimuth  of  the  reflecting  panel  facets  (the  predominant 
scene  component)  was  aligned  45°  and  135°  from  the  principal  plane 


Characteristic  Curves  (All) 


Figure  4.2-2.  (Note  that  the  dn  range  is  0  to  +510.) 

(containing  the  sun).  The  net  result  was  a  minimum  reflected  total 
radiance  among  the  five  scenes,  and  the  photo  CD  contrast  optimization 
algorithm  compensated  by  applying  higher  digital  numbers  to  the  scanned 
pixels.  The  upper  curves  in  Figure  4.2-1  plot  the  effect  of  this 
compensation.  Conversely,  Scenes  1,  2,  and  4  projected  more  reflected 
radiance  because  the  mean  panel  facet  azimuths  were  aligned  in  parallel 
with,  and  perpendicular  to,  the  principal  plane,  reflecting  more  sunlight  to 
the  sensor.  In  turn,  photo  CD  contrast  optimization  applied  lower  digital 
numbers  (dn)  to  these  images;  the  lower  curves  plot  the  effect  of  this 
compensation. 
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Figure  4.2-3  plots  the  average  characteristic  curve  for  the  two  darker 
scenes  (8  images)  as  the  upper  trace;  the  curve  for  the  three  brighter 
scenes  (12  images)  as  the  lower  trace;  and  the  total  collection  average 
(20  images)  as  the  middle  curve. 


Characteristic  Curves  (Averages) 


Figure  4.2-3. 

There  was  one  significant  side  benefit  of  the  contrast  optimization 
algorithm:  linear  stretching  of  the  log(exposure)  curve.  The  algorithm 
design  provides  linear  optimization  of  the  scanned  pixel  intensity  so  that 
digital  number  (dn)  approximates  actual  pixel  exposure/intensity. 

Figure  4.2-4  plots  log(digital  number)  versus  reflection  density  to 
demonstrate  the  effect  of  linear  optimization.  A  representative  film 
characteristic  curve  would  demonstrate  an  exponentially  increasing  region 
as  density  approaches  zero,  and  an  exponentially  decreasing  region  as 
density  aproaches  maximum  density.  Instead,  the  log  curve  assumes 
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near-linearity  across  the  range  of  measured  densities.  The  one  limitation 
of  this  plot  is  the  lack  of  measurement  for  higher  reflected  densities,  since 
the  lower  dn  region  is  not  mapped.  However,  the  trend  at  D  =  1.9  does  not 
indicate/predict  exponentially  decreasing  values  of  log(digital  number). 


Characteristic  Curves  (All  Images) 


Figure  4.2-4. 

A  least-squares  linear  fit  was  applied  to  all  the  uncalibrated  data  and  is 
plotted  on  Figure  4.2-4.  Since  normal  characteristic  curves  plot  density  as 
a  function  of  log(exposure),  the  mean  effective  system  gamma  (y )  is  the 
negative  reciprocal  of  the  fitted  line  slope  on  this  plot: 


AZ)  ^{reflection  density)  _  -1 

Alog^  A\og{digital  number)  -0.3306 


Equation  4-1 


In  lieu  of  more  comprehensive  sensitometric  data,  this  least-squares  fit 
line  was  used  as  the  common  baseline  for  calibrating  the  image  datasets, 
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i.e.,  the  pixel  digital  numbers  for  each  image  quad  were  interpolated  so 
that  the  least-squares  D-log(DN)  line  for  each  photo  pair  matched  this 
baseline.  Since 


[y  cal  )  -  bias^^, ]=D  =  ff  log(DN„„^^, )  -  bias^^^, ] 


Equation  4-2 


then 


ry  —  1  Q  uncal  log(  )-^bias^^,-bias^„^^i  )/y  ] 

cal 


Equation  4-3 


Image  analysis  of  the  radiometric  calibration 


The  calibrated  image  quads  were  used  to  create  the  Stokes  vector 
and  Stokes  derivative  images.  However,  since  the  sum  of  each  of  the  two 
orthogonally  polarized  image  pairs  in  the  quad  represents  an  estimate  of 
the  unpolarized  intensity  image  (I),  the  statistics  of  the  difference  between 
the  two  sums  can  provide  a  system  calibration  check  since,  if  the  system 
was  perfectly  calibrated,  the  difference  image  (D)  should  be  everywhere 
zero.  Reprising  Equation  3-6, 


D  =  r  - 1”  =  [|(0“)  +  1(90°)]  -  [1(45°)  +  1(135°)]  =  ~0 


Equation  4-4 


As  a  demonstration  of  calibration  performance,  post-calibration 
means  and  standard  deviations  were  calculated  for  66  x  66  pixel  regions 
within  the  imaged  neutral  density  card  visible  on  all  20  collections.  Figure 
4.2-5  plots  both  statistics  for  the  difference  (D)  subimages. 
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Figure  4.2-5.  (Note  that  the  full  dn  range  for  difference  images  is  -510  to  +510.) 

The  mean  difference  for  the  20  images  is  +3.4  dn  over  a  possible  range  of 
-510  to  +510  dn,  translating  to  a  bias  of  less  than  +0.7%.  The  mean  intra¬ 
image  (I’-l”)  standard  deviation  for  the  20  images  shows  a  fairly  constant 
value  of  about  6.6  dn,  while  the  mean  inter-image  standard  deviation  is 
about  one  half  that  amount,  3.3  dn.  Since  the  characteristic  curves  within 
each  image  quad  do  not  differ  as  much  as  they  do  between  independent 
image  collections,  this  difference  suggests  that  inter-image  calibration  was 
successful  at  reducing  variability  between  images  while  intra-image  (I’  vs. 
I”)  calibration  preserved  an  inherent  system  variation. 

As  a  second  demonstration,  post-calibration  means  and  standard 
deviations  were  calculated  for  the  same  10x10  pixel  subsets  from  the  10- 
step  wedge  visible  on  all  20  collections.  Figures  4.2-6  and  7  plot  these  two 
statistics  for  all  20  difference  (D)  subimages.  Figure  4.2-8  plots  the 
combined  statistics  for  the  20  difference  subimages. 


The  mean  differences  (Figure  4.2-6)  demonstrate  the  same  general 
positive  bias  as  was  previously  shown  in  the  neutral  density  card  data, 
with  digital  numbers  generally  ranging  between  -5  and  +10  dn  and  the 
combined  averages  (Figure  4.2-8)  showing  a  trend  from  +6  dn  near  D  =  0 
to  +2  dn  near  D  =  1.9. 

The  intra-image  (I’  - 1”)  standard  deviations  (Figure  4.2-7)  appear  to 
correlate  with  the  same  trend  for  intra-image  means:  the  combined 
deviation  decreasing  from  9-10  dn  near  D  =  0  to  5  dn  near  D  =1 .9.  The 
very  low  deviations  at  D  =0  correspond  to  intensity  values  that  were  close 
to  saturation  (dn  =  510)  for  those  image  subsets,  limiting  the  upper  range 
of  dn  variation. 
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Std  Dev  (Digital  Number) 


Calibrated  Step  Wedge  Standard  Deviation 
(All  Difference  Images) 


Figure  4.2-7. 


Calibrated  Step  Wedge  Statistics 
(Combined  Difference  Images) 


„  Average 


Reflection  Density 


Figure  4.2-8 


The  combined  inter-image  and  intra-image  standard  deviations  (Figure 
4.2-8)  compare  with  the  values  for  the  neutral  density  card:  the  intra¬ 
image  variation  is  consistently  (1 .5  to  2  times)  greater  than  the 
corresponding  inter-image  variation  at  a  given  value  of  reflection  density. 

Since  the  characteristic  curves  are  now  normalized,  the  residuals 
from  calibration  alone  will  sum  to  zero  bias.  One  explanation  for  the 
consistently  positive  bias  (+3  to  +6  dn)  between  the  calibrated  I’  image 
pair  and  the  I”  image  pair  is  clearly  evident  in  Figure  4.2-9. 


Calibrated  Step  Wedge  Average 
(Combined  Difference  Images  by  Solar  Elevation) 


»  1st  Sol  Elv 
__p._.._2nd  Sol  Elv 

. . 3rd  Sol  Elv 

■  . 4th  Sol  Elv 


Figure  4.2-9.  (Note  that  the  full  dn  range  is  -510  to  +510.) 


The  five  images  collected  at  each  of  the  four  solar  elevations  were 
combined  to  plot  their  mean  difference  over  reflection  density.  Since  each 
polarized  image  pair  originally  formed  one  35-mm  frame,  the  photo  CD 
contrast  optimization  algorithm  operated  on  the  pair  as  if  it  were  a  single 
image.  This  is  a  significant  problem  for  an  image  pair  that  has  high 
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contrast  in  one  member  and  low  contrast  in  the  other  member.  In 
particular,  this  condition  describes  the  0790°  image  pair,  which  maintains 
alignment  with  the  principal  plane  in  order  to  optimize  differential  contrast. 

For  image  quads  with  reduced  illumination,  e.g.,  the  4th  solar  elevation 
dataset,  the  overall  contrast  difference  within  the  two  pairs  would  be  small 
and  contrast  optimization  would  produce  very  similar  characteristic  curves. 
The  net  result  is  a  difference  (I’  - 1”)  image  which  approaches  zero,  as  is 
seen  in  Figure  4.2-9. 

However,  for  image  quads  with  increased  illumination,  e.g.,  the  1st  solar 
elevation  dataset,  the  overall  contrast  difference  would  be  larger  and 
contrast  optimization  would  produce  dissimilar  characteristic  curves. 
Although  a  normalized  characteristic  curve  is  applied  individually  to  each 
member  of  the  quad,  the  net  result  is  a  difference  image  which  preserves 
the  nonlinearities  of  the  differential  contrast  optimization  created  during 
photo  CD  processing. 

In  sum,  the  overall  effect  of  contrast  optimization  on  the  radiometric 
calibration  of  the  image  quads  is  a  balance.  However,  analysis  of  the 
calibration  data  presented  so  far  suggests  that  the  encountered  errors  are 
both  reasonable  and  explainable  if  not  removeable  or  reducible.  The  most 
significant  feature  of  the  calibration  is  the  presence  of  a  systematic  bias 
between  the  estimates  of  I’  and  I”  that  varies  between  +0.5%  and  +1.0%. 

Polarizer  reference  alignment  analysis 

Minor  rotations  of  the  exterior  orientation  parameter  kappa  (k) 
were  visible  in  the  imagery.  (Refer  to  Figure  4.1-1 .)  The  mounting  frame 
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for  the  polarimeter  apparently  underwent  a  small  amount  of  operator- 
induced  twisting  each  time  the  film  advance  lever  was  manipulated.  The 
amount  of  the  rotation  is  plotted  in  Figure  4.2-10. 


Polarizer  Reference  Misalignment  Measurement 


Figure  4.2-10.  (Note  that  the  full  rotation  range  is  -180°  to  +180°.) 

The  procedure  for  measuring  this  rotation  was  to  measure  the  pixel 
coordinates  of  the  two  endpoints  of  the  line  segment  described  by  the 
abutment  of  the  two  test  panels.  Since  the  panels  always  maintained 
relative  alignment  with  the  reference  plane  in  increments  of  45°,  the 
rotational  alignment  error  of  the  polarizer  relative  to  the  reference  plane 
could  be  measured  for  each  image: 

Rotation  Error  =  tan''|  —  ]  +  «(45"),  n  =  (0,1, 2, 3} 

\AxJ 


Equation  4-5 
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The  results  show  that  the  polarizer  was  positively  rotated  over  a 
range  between  +1 .4°  and  +1 .6°  during  the  first  solar  elevation  collection 
and  was  significantly  increased  to  a  range  between  +5.5°  and  +7.4° 
during  subsequent  collections. 

Analysis  of  the  mirror  image  subsets  will  demonstrate  both  the  effect  of 
this  rotation  error  on  the  visualization  of  the  reflected  polarized  skydome 
and  the  post-processing  capability  to  compensate  the  Stokes  imagery  for 
this  alignment  error. 

Solar  elevation  measurement  analysis 

The  length  of  the  shadow  formed  by  a  small  vertical  post  mounted  on  the 
platform  was  used  to  calculate  the  solar  elevation  angle.  The  results  of 
this  measurement  are  listed  in  Table  4-2. 


Table  4-2 

Solar  Elevation  Secondary  Measurement  Results 


Dataset  Nbr 

Solar  Elevation 

Std  Dev 

1st 

+26.9° 

0.2° 

2nd 

+23.0° 

0.06° 

3rd 

+16.3° 

0.04° 

4th 

+7.9° 

0.04° 

These  values  are  plotted  in  Figure  4.2-11  against  the  calculated  solar 
elevation  angle  versus  local  time.  Sunrise  (0°)  was  at  0645  EST  and 
sunset  (0°)  was  at  1650  EST. 


These  measured  solar  elevation  values  were  used  as  inputs  within  the 
imaging  simulation. 


97 


Histogram  analysis  of  the  intensity  (H  full  images 


Figure  4.2-12  plots  all  20  histograms  for  the  calibrated  intensity  (I) 
full  images.  This  plot  and  subsequent  plots  demonstrate  the  high 
variability  of  intensity  (I)  relative  to  the  other  two  Stokes  parameters  (P 
and  T).  This  variability  is  directly  attributable  to  the  effects  of  imaging 
geometry,  as  will  be  shown  in  Figures  4.2-13  and  14. 

There  are  two  general  features  to  note  in  these  figures: 

First,  the  two  test  panels  are  highly  separable  in  the  I  histograms.  The 
aggregate  of  the  dark  panel  pixels  have  dn  values  centered  in  the 
neighborhood  of  dn  =  200  and  the  light  panel  pixels  are  aggregated 
around  dn  =  400.  The  dark  panel  pixel  distributions  have  about  one  half 
the  variance  of  the  dark  panel  distributions.  (The  full  statistical  comparison 
of  the  panels  will  be  treated  during  sub-image  analysis  which  follows.) 

Second,  the  combined  effects  of  film  development  and  photo  CD  contrast 
optimization  are  apparent  at  the  extrema.  Almost  no  darker  pixels  are 
mapped  to  the  lower  dn  values  between  0  and  40.  This  apparent  failure  of 
contrast  optimization  is  instead  due  to  the  presence  of  unexposed  film 
area  within  the  processed  35-mm  frame.  Because  of  the  smaller  exposure 
area  in  the  Nishika  camera,  approximately  9%  of  a  standard  35-mm  frame 
is  unexposed  and  the  resulting  processed  pixels  have  dn  values  of  zero. 
These  dark  pixels  are  considered  in  the  contrast  optimization,  effectively 
adding  a  positive  bias  of  20  dn  to  the  darkest  exposed  pixels  (40  dn  in  the 
I  image).  The  3-5%  of  saturated  pixels  (dn  values  approaching  510)  are 
the  result  of  actual  film  saturation,  contrast  optimization,  or  both.  These 
saturated  pixels  have  a  significant  effect  in  the  histogram  analysis  of  the  T 
(polarization  ellipse  orientation  angle)  images  which  follows. 
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Figure  4.2-13  illustrates  the  variation  of  intensity  relative  to  solar 
elevation  angle.  In  this  case,  reflected  intensity  generally  increases  as  the 
solar  elevation  angle  decreases,  down  to  a  minimum  that  appears  to  be 
near  the  4th  solar  elevation.  The  key  effect  is  due  to  the  net  Fresnel 
reflectivity  of  the  panel  facet  slope  distribution.  The  small  but  increasing 
distribution  neardn  =  100  are  pixels  in  shadow. 

Figure  4.2-14  illustrates  the  variation  of  intensity  relative  to  scene 
geometry.  In  scenes  1,  2,  and  4,  the  mean  azimuth  angle  of  the  panel 
facet  distribution  is  aligned  with  the  principal  plane,  favoring  reflection  of 
direct  sunlight.  In  scenes  3  and  5,  the  mean  azimuth  angle  is  rotated  -45° 
and  +45°,  respectively,  from  the  principal  plane.  The  main  benefit  of  using 
textured  panels  is  that  the  reflection  geometry  of  the  scene  is  both  a 
deterministic  and  controllable  feature  of  the  imaging  experiment. 
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Histogram  Analysis  of  the  %  polarization  CP)  full  images 


Figure  4.2-15  plots  all  20  histograms  for  the  calibrated  % 
polarization  (P)  full  images.  This  plot  and  subsequent  plots  demonstrate 
both  the  relatively  low  variability  of  %  polarization  (P)  compared  with 
intensity  (I)  and  its  relative  independence  of  imaging  geometry. 


The  most  significant  feature  of  the  P  histograms  is  the  relative 
absence  of  pixels  with  polarizations  approaching  100%:  on  average,  less 
than  0.5%  of  all  the  pixels  in  the  20  images  have  polarizations  above 
65%.  The  mode  value  averages  1 1 .8%  with  a  standard  deviation  of  0.9%. 
In  the  absence  of  an  independent  system  for  calibrating  the  polarization 
measurements  to  an  absolute  scale,  the  combined  results  indicate  that 
the  measurements  are  statistically  consistent  over  a  relative  scale.  The 
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P  Histogram 

(Combined  Scene  Full  Images) 
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Figure  4.2-16.  (Note  the  reduced  scale  for  %  polarization.) 


expectation  for  non-ideal  polarimeter  performance  is  to  observe  some 
amount  of  systemic  depolarization,  primarily  due  to  partial  transmission  of 
the  orthogonal  radiation  component  through  a  non-ideal  linear  polarizer. 
However,  because  there  are  four  intensity  measurements  of  radiation 
passing  through  the  same  non-ideal  polarizing  material,  the  net  observed 
depolarization  should  be  proportional  for  each  Stokes  component.  In  the 
calculation  of  polarization,  this  depolarization  factor  appears  in  both  the 
numerator  and  denominator  (Equation  3-4)  and  thus  cancels  out. 


Figures  4.2-16  and  17  do  not  demonstrate  any  observable  polarization 
trend  with  respect  to  solar  azimuth  or  scene  geometry.  A  more  detailed 
analysis  of  individual  features  is  provided  in  Chapter  4.5. 
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P  Histogram 

(Combined  Soiar  Elevation  Full  Images) 
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Figure  4.2-17.  (Note  the  reduced  scale  for  %  polarization.) 


Histogram  Analysis  of  the  orientation  angle  (T)  full  images 


Figure  4.2-18  plots  all  20  histograms  for  the  calibrated  orientation 
angle  (T)  full  images.  This  plot  and  subsequent  plots  demonstrate  the 
sensitivity  of  T  to  the  presence  of  saturated  pixel  intensities  as  well  as  the 
rotational  variation  of  kappa  (k  ).  (Refer  to  Figure  4.1-1.) 


The  significant  spikes  at  22.5°  intervals  are  mathematical  artifacts 
created  by  the  presence  of  one  or  both  image  pairs  containing  pixels  with 
saturated  intensity  values,  i.e.,  dn  =  255.  Referring  back  to  Equation  3-5, 
T  is  the  half-angle  of  the  arctangent  of  the  ratio  of  U  and  Q.  The  high  pixel 
pixel  count  spikes  occur  at  2xT  values  of  +/-  0°,  45°,  90°,  135°,  and  180°. 


104 


Orientation  Angle  (degrees) 


Figure  4.2-18. 


The  tangent  of  these  angles  corresponds  to  the  following  possible 
combinations  of  U  and  Q: 


Table  4-3 

U  &  Q  Combinations 


Angle 

Tan(Angle) 

Possible  U  &  Q  values 

0“ 

0 

U  =  0,  Q  =  {0..+510} 

+45° 

+1 

U  =  Q  =  {+1..+510} 

+90° 

+  00 

U  =  {+1..+510},  Q  =  0 

+135° 

-1 

-U  =  Q,  U  =  |Q|  =  {+1..+510} 

±180° 

0 

U  =  0,  Q  =  {-1..-510} 

-135° 

+1 

U  =  Q  =  {-1..-510} 

-90 

-00 

U  =  {-1.-510},  Q  =  0 

-45° 

-1 

U  =  -Q,U  =  |Q|  =  {-1.-510} 

Figure  4.2-19. 


TC  Histogram 

(Combined  Solar  Elevation  2+3+4  Full  Images) 
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Figure  4.2-20. 
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The  largest  spikes  occur  where  either  U  or  Q  are  0;  the  single 
largest  spike  occurs  when  U  and  Q  are  both  0.  The  less  prominent  spikes 
occur  when  U  and  Q  fluctuate  between  values  of-1  and  +1.  This 
incidence  is  observed  for  near  and  fully  saturated  pixel  intensities  within 
each  polarized  pair  of  an  image  quad.  The  same  situation  would  exist 
near  zero  intensity.  However,  this  situation  was  not  observed  within  this 
experiment  because  of  the  large  positive  dn  bias  identified  in  the  I  image 
histogram  analysis. 

If  these  spikes  are  disregarded,  the  remaining  distribution  of  T 
values  provides  a  unique  indicator  for  the  azimuth  location  of  the  principal 
plane  relative  to  the  polarimeter  reference  plane.  Figure  4.2-19  illustrates 
the  combined  average  histogram  of  the  first  solar  elevation  collection, 
where  the  polarimeter  rotation  bias  is  approximately  +1.5°;  and  Figure  4.2- 
20  combines  the  remaining  collections,  where  the  bias  is  in  the  vicinity  of 
+6.5°.  The  distribution  of  T  values  approach  a  maximum  pixel  count  at  the 
azimuth  location  of  the  principal  plane.  This  distribution  correlates  with  the 
attenuated  reflection  of  skylight  from  a  surface  containing  reflecting  facets 
that  assume  a  near-uniform  distribution  of  azimuths  and  a  near-Gaussian 
distribution  of  slopes,  for  example,  a  waterwave  surface.  In  fact,  this  same 
distribution  is  visible  for  the  water  surface  T  image  histogram  analyzed  in 
Chapter  4.4. 

In  sum,  the  existence  of  truncated  values  at  either  extremum  of  the 
intensity  scale  will  produce  mathematical  artifacts  in  the  generation  of  the 
T  parameter  image,  creating  biases  in  the  estimate  of  mean  and  variance. 
In  the  absence  of  these  artifacts,  the  orientation  angle  of  maximum  pixel 
count  can  serve  as  an  indicator  of  the  principal  plane  azimuth  in  imaged 
scenes  with  near-homogeneous  surface  geometry  statistics. 
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4.3  Image  Simulation  Results 


The  intent  of  this  qualitative  evaluation  is  compare  the  relative 
spatial  distributions  of  clear  sky  Stokes  parameter  measurements  from 
actual  images  of  a  reflecting  dome  mirror  with  the  predictions  of  an 
imaging  simulation  that  models  the  collection  geometry.  With  reference  to 
Chapter  3,  no  attempt  was  made  to  calibrate  either  the  model  or  the 
imaging  system  for  specific  spectral  response  or  absolute  radiometry.  The 
benefit  of  this  qualitative  approach  is  to  provide  a  reasonableness  check 
on  the  spatial  performance  of  the  system  and  model,  consistent  with  the 
objectives  of  this  study.  Chapter  5  provides  several  suggestions  for  further 
research  which  would  seek  to  both  validate  the  model  and  calibrate  the 
system. 

With  reference  to  Figures  4.3-1  through  4,  the  one  signficant 
deviation  is  the  existence  of  image  saturation  in  the  vicinity  of  the  sun  and 
the  lack  of  same  in  the  simulation. 

With  reference  to  Figures  4.3-5  through  8,  the  image  and 
simulation  both  demonstrate  the  shifting  of  the  peak  polarization  band  that 
is  orthogonal  to  the  solar  point  (i.e.,  the  anti-solar  point).  Again,  saturation 
is  evident  in  the  vicinity  of  the  sun.  However,  the  sun,  being  an  fully 
unpolarized  source,  would  otherwise  be  represented  as  dark. 

With  reference  to  Figures  4.3-9  through12,  the  darkest  pixels 
correspond  to  T  values  approaching  +90°;  the  lightest  pixels  correspond 
to  T  values  approaching  -90°.  The  most  significant  deviation  with  these 
figures  is  due,  unfortunately,  to  the  color  and  gray  scale  representations  of 
two  different  software  programs  (ENVI™  and  MathCAD™);  the  color 
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renderings  provide  better  spatial  correlation.  However,  the  model  is 
relatively  insensitive  to  the  rotational  bias  that  was  also  simulated  to 
coincide  with  the  later  three  solar  elevation  collections. 

With  reference  to  Figures  4.3-13  throughi 6,  the  abs(T)  images 
clearly  demonstrate  the  sensitivity  of  the  imaging  system  to  rotation  error. 
The  first  abs(T)  image  provides  an  example  of  the  bi-lateral  symmetry  that 
is  expected  for  abs(T)  values  on  hemispheres  bisected  by  the  principal 
plane. 

The  predominant  feature  in  the  image  that  demonstrates  this  symmetry  is 
the  dark  (near-0°)  areas  opposing  the  hemisphere  containing  the  sun.  The 
presence  of  the  polarimeter  supports  in  the  scene  provides  a  source  of 
reference  for  the  detection  of  symmetry  (or  lack  thereof). 

With  reference  to  Figures  4.3-17  through  20,  both  the  images  and 
their  simulations  have  undergone  a  mathematical  rotation  with  the 
appropriate  bias  applied,  via  atan(tan(T-bias)).  The  signficant  feature  of 
this  rotation  is  the  approximate  restoration  of  bi-lateral  symmetry  for  all 
four  images. 

The  dark  (near-0°)  areas  all  demonstrate  relative  symmetry  with  respect  to 
the  polarimeter  supports,  as  would  be  expected  for  polarization  due  to 
pure  Rayleigh  scattering  under  clear  sky  conditions. 
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Figures  4.3-1  through  4  (a)  and  (b).  I  parameter  images  and  their  simulations 
for  solar  elevations  26.9°,  23,0°,  16.3°,  and  7.9°,  respectiveiy. 
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Figures  4.3-5  through  8  (a)  and  (b).  P  parameter  images  and  their  simulations 
for  solar  elevations  26.9°,  23.0°,  16.3°,  and  7.9°,  respectively. 
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Figures  4.3-9  through  12  (a)  and  (b).  T  parameter  images  and  their  simulations 
for  solar  elevations  26.9°,  23.0°,  16.3°,  and  7.9°,  respectively. 


ATMAP 


ATMAP 


Figures  4.3-13  through  16  (a)  and  (b).  Absolute  value  of  T  parameter  images  and  their 
simulations  for  solar  elevations  elevations  26.9°,  23.0°,  16.3°,  and  7.9°,  respectively.  Note 
the  asymmetry  of  the  spatial  distributions  of  T  compared  with  Figures  4.3-17  through  20. 


Figures  4.3-17  through  20  (a)  and  (b).  Absolute  value  of  T  parameter  images  and  their 
simulations  after  having  a  rotational  bias  of  -1 .4°,  -5.9°,  -5.6°,  and  -6.2°  (respectively) 
applied  to  each  pixel  value,  corresponding  to  the  rotation  error  of  the  polarimeter  relative 
to  the  principal  plane.  Note  the  symmetry  of  the  spatial  distributions  of  T  compared  with 
Figures  4.3-13  through  16. 


Although  unmodeled  by  the  simulation,  the  mirror  sub-scene 
difference  images  are  presented  as  Figures  4.3-21  through  24.  Within  the 
context  of  potential  spatial  system  calibration,  these  images  provide 
insight  into  the  underlying  accuracy  of  the  intensity  measurements  that 
were  used  to  calculate  the  Stokes  parameter  images  and  their  derivative 
images. 

Although  the  in-field  test  scenes  were  relatively  flat  compared  with 
the  scale  of  focal  distance,  one  observable  instance  of  misregistration  is 
evident  in  the  D  images:  the  largest  deviations  in  the  surface  plots 
correspond  to  misregistrations  of  the  10-foot  polarimeter  supports  from 
parallax  effects  created  by  the  finite  separation  of  the  four  apertures  (55 
mm  between  the  two  most  distant  lens  centers).  The  effect  is  especially 
pronounced  at  the  outer  edges  of  the  dome  mirror  where  secondary 
reflection  of  the  polarimeter  supports  enhanced  the  dispersion  of  their 
reflected  images. 
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X  =  +7.8051 
a  =  20.3882 


X  =  +3.7067 
a  =  16.2365 


4.4  Aerial  calibration  results 


General  collection  results 


Twenty  photo  quads  (40  35-mm  frames)  were  collected  between 
1312  and  1548  EST  on  03  Nov  1994.  A  summary  of  the  collection  and 
processing  results  is  listed  in  Table  4-4. 


Table  4-4 

Collection  and  Processing  Results  Summary 


Time 

EST 

Ht 

ft 

Az 

deg 

Quad 

Nbr 

Frame 

Nbrs 

Remarks 

1312 

1000 

045 

1 

16/- 

water  /*partially  processed 

1313 

1000 

045 

2 

-/1 7 

water  /*partially  processed 

1314 

500 

090 

3 

18/- 

water  /*partially  processed 

1314 

500 

090 

4 

water  /  unprocessed 

1315 

150 

270 

5 

water  /  unprocessed 

1315 

150 

270 

6 

water  /  unprocessed 

1319 

1000 

unk 

7 

19/20 

target  of  opportunity 

1320 

1000 

unk 

8 

21/22 

target  of  opportunity 

1323 

1000 

unk 

9 

23/24 

target  of  opportunity 

1327 

1000 

unk 

10 

25/26 

target  of  opportunity 

1530 

1000 

045 

11 

water  /  unprocessed 

1530 

1000 

045 

12 

water  /  unprocessed 

1531 

500 

090 

13 

27/28 

water  /  processsed 

1532 

500 

090 

14 

water  /  unprocessed 

1533 

150 

270 

15 

water  /  unprocessed 

1533 

150 

270 

16 

water  /  unprocessed 

1536 

1000 

unk 

17 

30/31 

target  of  opportunity 

1538 

1000 

unk 

18 

32/33 

target  of  opportunity 

1545 

1000 

unk 

19 

34/35 

target  of  opportunity 

1548 

1000 

unk 

20 

36/37 

target  of  opportunity 

*  Photo  CD  processing  accepted  one  quad  frame  but  rejected  the  other. 


As  previously  mentioned  in  Chapter  4.1,  the  aerial  collection 
experienced  two  failures  as  a  direct  result  of  photo  CD  processing. 
Because  the  contrast  optimization  algorithm  was  not  disabled  during 
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photo  CD  processing,  the  sensitometric  step  wedge  was  both  distorted 
and  incomplete,  resulting  in  the  loss  of  sensitometric  calibration.  Also,  all 
but  one  of  the  12  Hinckley  reservoir  collections  were  partially  or  totally 
rejected  by  the  optimization  algorithm  for  failure  to  meet  a  minimum 
threshhold  for  input  contrast.  The  one  quad  that  was  fully  processed 
exceeded  the  threshhold  only  because  it  contained  a  small  amount  of 
brighter  land  area  with  the  otherwise  dark  waterwave  surface  in  the 
scene,  providing  a  range  of  contrast  that  was  acceptable  to  the 
processing  algorithm. 

Due  to  the  overall  failure  of  the  overwater  collections,  the  one  fully 
processed  quad  was  used  as  a  demonstration  of  feasibility  for  further 
aerial  polarimetric  imaging  activities  based  upon  its  general  agreement 
with  the  results  of  the  in-field  experiment.  The  most  significant  loss  to  this 
demonstration  is  the  absence  of  sensitometric  calibration. 

However,  there  are  several  unique  aspects  of  the  one  aerial 
collection  that  can  still  be  evaluated.  In  particular,  the  relative 
homogeneity  of  the  waterwave  surface  provides  a  basis  for  comparison  of 
the  spatial  information  content  between  contiguous  sub-scenes.  In  this 
feasibility  demonstration,  one  large  (512  x  512  pixel)  sub-image  of  the 
watenwave  surface  was  subdivided  into  four  256  x  256  pixel  contiguous 
patches.  The  sub-image  statistics  could  then  be  compared  under  the 
assumption  that  the  four  sub-scenes  represent  four  correlated  realizations 
of  the  same  underlying  process,  i.e.,  wind-driven  waterwave  generation. 

Figures  4.4-1  through  4  depict  the  512  x  512  pixel  Stokes  derivative  sub¬ 
images.  Again  I  is  the  unpolarized  intensity  image,  P  is  the  %  polarization 
image,  T  is  the  image  depicting  the  orientation  angle  of  the  polarization 
ellipse,  and  D  is  the  difference  image,  representing  I’  - 1”. 
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Patch  B  I  Patch  D 

(256x256  1(256x256 
pixels)  I  pixels) 


Figure  4.4-5. 


Histogram  analysis  of  the  difference  (D)  image 


Figure  4.4-6  plots  the  four  histograms  for  the  uncalibrated  difference  (D) 
image  patches. 
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Figure  4.4-6. 


In  lieu  of  sensitometric  calibration,  this  plot  provides  an  indication  of  the 
magnitude  of  error  that  exists  in  the  calculation  of  the  Stokes  parameters 
and  Stokes  derivatives. 

The  underlying  statistics  of  the  photo  quad  for  the  512  x  512  sub-image 
and  four  256  x  256  patches  are  presented  in  Table  4-5. 
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Parameter 

Min 

Table  4-5. 
Statistics  Report 
Max 

Mean 

Std  Dev 

1  (0°)  (Subimage) 

81.0 

255.0 

179.4123 

38.7929 

1  (90°)  (Subimage) 

5.0 

246.0 

88.3622 

25.6739 

1  (45°)  (Subimage) 

22.0 

248.0 

130.0861 

35.9600 

1  (135°)  (Subimage)  83.0 

255.0 

186.7682 

38.5032 

D  (Subimage) 

-173.0 

124.0 

-49.0798 

26.3835 

D  (Patch  A) 

-147.0 

55.0 

-41.6093 

22.7740 

D  (Patch  B) 

-173.0 

25.0 

-69.1875 

21.8237 

D  (Patch  C) 

-148.0 

124.0 

-30.9919 

22.8967 

D  (Patch  D) 

-146.0 

34.0 

-54.5306 

21.1701 

The  most  significant  feature  of  the  uncalibrated  difference  image 
patches  is  the  strong  negative  bias,  between  -30  and  -70  digital  numbers 
(dn).  These  differences  are  systematic  errors  created  by  the  absence  of 
sensitometric  calibration  and  the  presence  of  photo  CD  contrast 
optimization.  As  will  be  seen  in  the  T  image  data,  the  reference  azimuth  is 
in  the  vicinity  of  the  principal  plane.  The  high  contrast  between  the  I  (0°) 
and  I  (90°)  image  patches  bear  out  this  fact. 


However,  the  dn  values  for  the  I  (45°)  and  I  (135°)  sub-images  appear  to 
be  artificially  elevated  compared  with  the  other  cross-polarized  pair  and 
their  expected  conformance  with  the  law  of  Malus:  the  mean  for  the  I  (45°) 
sub-image  has  an  expected  dn  value  of  109  and  the  mean  for  the  I  (135°) 
has  an  expected  dn  value  of  158.  But  because  this  polarized  pair  has 
overall  lower  exposure  intensity,  contrast  optimization  applied  a  significant 
positive  dn  bias  of  +16%  to  the  pair  within  the  image  frame  and  created  a 
difference  (□  =  !’- 1”)  of  49  dn  between  the  two  cross-polarized  pairs. 


Another  analytic  feature  of  this  uncalibrated  data  is  the  variation  of 
the  patch  differences.  Figure  4.4-4  illustrates  the  spatial  distribution  of 
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differences  values  across  the  four  patches  in  the  sub-image.  Because  the 
patch  differences  are  all  the  product  of  a  single  contrast  optimization 
application,  they  are  created  by  differential  contrast  optimization  applied  to 
the  pair  within  each  frame.  The  effect  of  this  non-linear  optimization  also 
modifies  the  calculation  of  Q  and  U,  creating  subsequent  errors  in  the 
calculation  of  P  and  T. 

In  sum,  the  uncalibrated  image  has  a  systematic  error  of -9.6%. 
While  this  bias  is  over  ten  times  the  magnitude  of  the  calibrated  images, 
the  only  alternative  is  to  reject  this  one  remaining  image  from  the  aerial 
collection  for  lack  of  collateral  sensitometric  calibration  data.  A  future 
option  for  consideration  is  to  utilize  known  characteristic  curves  for  similar 
histogram  content  and  execute  an  iterative  least  squares  correction  of  the 
photo  quads  on  a  pixel  by  pixel  basis  until  the  resulting  difference  image 
achieves  an  acceptable  minimum.  This  type  of  correction  is  currently 
outside  the  range  of  ENVI  processing  capabilities;  however,  the  merits  of 
this  option  will  be  addressed  in  Chapter  5. 

Histogram  analysis  of  the  unpolarized  intensity  (H  image 

Figure  4.4-7  plots  the  four  histograms  for  the  uncalibrated  intensity  (I) 
image  patches.  In  lieu  of  sensitometric  calibration,  this  plot  provides  only 
an  indication  of  relative  intensities,  with  a  known  systematic  error  of 
between  +15  and  +35  dn  in  the  I  parameter  calculation. 

The  underlying  statistics  of  the  uncalibrated  512  x  512  sub-image  and  four 
256  X  256  patches  are  presented  in  Table  4-6. 


122 


I  Histogram  (uncalibrated) 
(4  contiguous  patches) 
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Digital  Number 


Figure  4.4-7. 


Parameter 

Min  Max 

Table  4-6. 
Statistics  Report 
Mean 

Exp.  Mean 

Std  Dev 

I  (Subimage)  139.0  457.0 

292.3145 

267.7746 

64.2228 

1  (Patch  A) 

143.0  335.0 

233.2161 

212.4115 

27.0159 

1  (Patch  B) 

210.5  445.5 

329.4756 

294.8819 

39.4741 

1  (Patch  C) 

139.0  368.0 

246.4961 

231.0002 

30.8259 

1  (Patch  D) 

248.0  457.5 

360.0700 

332.8047 

41.1861 

The  figure  and  table  again  illustrate  the  large  variation  of  I  relative 
to  P  and  T;  in  this  example,  the  inter-patch  standard  deviation  is  about  62 
dn,  comparable  to  the  full  subimage  standard  deviation.  If  the  expected 
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means  are  used  as  estimates  of  the  calibrated  I  means,  the  inter-patch 
standard  deviation  decreases  by  only  6  dn,  to  56  dn. 

Figure  4.4-1  illustrates  the  presence  of  strong  solar  reflection  in  Patches  B 
and  D,  providing  strong  contrasts  to  Patches  A  and  C  and  also  providing 
an  explanation  for  the  non-systematic  inter-patch  variation.  While  the 
strong  reflection  of  sunlight  is  an  inadvertant  byproduct  of  yet  another 
camera  orientation  error,  deviations  from  a  non-level  platform  (errors  in 
omega  or  phi,  refer  to  Figure  4.1-1),  it  provides  a  more  dramatic 
comparison  with  the  relative  stability  of  the  P  and  T  images  and  their 
spatial  spectra. 

Histogram  analysis  of  the  %  polarization  (PI  image 

Figure  4.4-8  plots  the  four  histograms  for  the  uncalibrated  %  polarization 
(P)  image  patches.  In  lieu  of  sensitometric  calibration,  this  plot  provides 
only  an  indication  of  approximate  polarizations,  with  an  unknown 
systematic  error  in  the  P  parameter  calculation  due  to  the  non-linearity  of 
differential  contrast  optimization. 

The  underlying  statistics  of  the  uncalibrated  512  x  512  sub-image  and  four 
256  X  256  patches  are  presented  in  Table  4-7. 

The  figure  and  table  illustrate  the  small  variation  of  P  relative  to  1;  in 
this  example,  the  inter-patch  standard  deviation  is  2.43  %,  less  than  one 
third  of  the  subimage  standard  deviation.  Also,  the  intra-patch  standard 
deviations  are  comparable  with  the  subimage  standard  deviation. 
However,  the  P  data  exhibits  the  same  spatial  trend  visible  in  the  mirror- 
reflected  skydome  image,  consistent  with  the  observation  that  the  water 
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surface  represents  a  quasi-Gaussian  distribution  of  reflecting  surface 
facets  that  reflect  a  smaller  solid  angle  portion  of  the  skydome  relative  to 
the  full  coverage  of  the  dome  mirror. 


Figure  4.4-2  illustrates  the  spatial  distribution  of  P  with  the  presence  of 
lower  values  in  the  vicinity  of  the  solarization  area  in  Patches  B  and  D  and 
higher  values  in  Patches  A  and  C  away  from  the  area  of  solar  reflection. 


125 


Consistent  with  the  skydome  reflection  results  and  the  simulation  model 
results  of  Chapter  4.3,  reflected  P  values  are  minimum  at  the  point  of 
direct  sun  reflection  and  maximum  in  the  area  of  antisolar  point  reflection. 

Histogram  analysis  of  the  orientation  angle  (T)  image 

Figure  4.4-9  plots  the  four  histograms  for  the  uncalibrated  orientation 
angle  (T)  image  patches.  In  lieu  of  sensitometric  calibration,  this  plot 
provides  only  an  indication  of  approximate  angles,  with  an  unknown 
systematic  error  in  the  T  parameter  calculation  again,  as  with  P,  due  to 
the  non-linearity  of  differential  contrast  optimization. 


T  Histogram  (uncalibrated) 
(4  contiguous  patches) 
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Figure  4.4-9. 


The  underlying  statistics  of  the  uncalibrated  512  x  512  sub-image  and  four 
256  X  256  patches  are  presented  in  Table  4-8. 
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Parameter 

Min 

Table  4-8. 
Statistics  Report 
Max 

Mean 

Std  Dev 

T  (Subimage) 

-82.62 

65.99 

-16.06 

6.00 

T  (Patch  A) 

-82.62 

65.99 

-15.12 

5.84 

T  (Patch  B) 

-32.97 

16.94 

-11.33 

4.46 

T  (Patch  C) 

-61.43 

18.31 

-21.02 

5.02 

T  (Patch  D) 

-39.06 

10.28 

-16.77 

4.06 

The  figure  and  table  illustrate,  as  with  P,  the  small  variation  of  T 
relative  to  1;  in  this  example,  the  inter-patch  standard  deviation  is  4.02°, 
less  than  two  thirds  the  subimage  standard  deviation.  Also,  the  intra-patch 
standard  deviations  are  again  comparable  with  the  subimage  standard 
deviation.  However,  the  T  data  exhibits  the  same  spatial  trend  visible  in 
the  mirror-reflected  skydome  image,  consistent  with  the  observation  that 
the  water  surface  represents  a  quasi-Gaussian  distribution  of  reflecting 
surface  facets  that  reflect  a  smaller  solid  angle  portion  of  the  skydome. 

Figure  4.4-3  illustrates  the  spatial  distribution  of  T.  The  flight  azimuth  and 
hence  polarizer  reference  for  this  collection  was  approximately  90°  (due 
East).  The  apparent  position  of  the  water-reflected  sun  in  the  image  is 
approximately  WSW.  Therefore,  the  polarizer  reference  plane  is  oriented 
from  the  principal  plane  by  the  approximate  azimuth  angle  of  -16°  as 
calculated  from  the  mean  of  the  T  subimage.  This  value  compares  closely 
with  the  T  calculation  using  the  photo  quad  means  in  Table  4-5: 


T  -  0.5  tan  '  =  =  0.5  tan  ' 

(130.0861-186.7682) 

=  -15.95° 

Equation  4-5 

Q 

(179.4123- 88.3622)  _ 

Patch  D  is  closest  to  the  location  of  the  solar  reflection  point  and,  as  a 
result,  contains  the  more  even  distribution  of  positive  and  negative  T 
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values:  the  principal  plane  roughly  bisects  the  patch  from  the  bottom  right 
corner  of  the  patch  and  runs  through  the  top  right  corner  of  Patch  A.  T 
values  in  Patch  B  are  more  positive  than  the  mean,  consistent  with  water 
surface  facets  reflecting  the  polarized  skydome  from  one  side  of  the 
principal  plane;  conversely,  values  in  Patch  C  are  more  negative  than  the 
mean,  consistent  with  water  surface  facets  reflecting  the  polarized 
skydome  from  the  other  side  of  the  principal  plane.  Patch  A  is  slightly  less 
positive  than  Patch  B  since  the  principal  plane  traverses  it. 

Consistent  with  the  polarizer  orientation  error  correction  applied  in 
Chapter  4.3,  the  -16°  bias  can  also  be  removed  from  this  image  through 
the  same  process.  Figure  4.4-10  and  Table  4-9  show  the  effects  of  this 
mathematical  rotation  of  the  orientation  angle  of  the  polarization  ellipse. 

In  particular,  the  means  and  distributions  of  Patches  D  and  A  are  nearly 
centered  on  zero  while  the  means  and  distributions  of  Patches  B  and  C 
are  biased  positive  and  negative  by  roughly  equal  amounts,  respectively. 

One  additional  observation  in  Figures  4.4-9  and  10  is  the  absence 
of  artifact  spikes  in  the  T  histograms:  saturated  dn  values  do  exist  in 
small  numbers  but,  from  Table  4-5,  it  is  apparent  that  they  do  not  occur 
together  in  either  polarized  pair.  Therefore,  the  values  of  Q  or  U  that  are 
zero  and  near-zero  occur  in  proportion  with  the  rest  of  the  histogram 
distribution. 
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T  Histogram  (corrected) 
(4  contiguous  patches) 
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Figure  4.4-10. 


Parameter 

T  (Subimage) 
T  (Patch  A) 

T  (Patch  B) 

T  (Patch  C) 

T  (Patch  D) 


Table  4-8. 
Statistics  Report 


Min 

Max 

Mean 

Std  Dev 

-65.86 

82.76 

+0.71 

6.00 

-65.86 

82.76 

+1.65 

5.84 

-16.20 

33.71 

+5.44 

4.46 

-44.66 

35.08 

-4.25 

5.02 

-22.29 

27.05 

+0.00 

4.06 

Figure  4.4-12.  Abs(T)  image. 
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Figure  4.4-11.  Abs(T  +  bias) 


As  in  Chapter  4.3,  Figures  4.4-11  and  12  illustrate  the  absolute 
value  of  T  before  and  after  the  addition  of  the  bias  correction.  Again,  the 
dark  pixels  correspond  to  surface  facets  with  azimuths  oriented  so  that 
light  is  reflected  from  the  principal  plane  where  T  is  near  0°;  and  lighter 
pixels  correspond  to  surface  facets  with  azimuths  increasingly  oriented 
away  from  the  principal  plane.  The  rotated  image  of  Figure  4.4-12 
correlates  favorably  with  the  sun  glint  pattern  in  Figure  4.4-1 ,  providing  a 
secondary  indication  of  the  reflected  principal  plane. 

In  the  absence  of  multiple  and  calibrated  waterwave  images,  the 
remainder  of  the  analysis  of  the  one  uncalibrated  image  falls  within  the 
description  of  a  feasibility  demonstration.  In  Chapter  4.5,  the  power 
spectral  densities  (PSDs)  of  the  four  patches  were  evaluated  for  their 
stability  and  compared  with  the  results  of  the  in-field  experiment.  While 
the  in-field  PSD  stability  analysis  evaluates  the  imaging  geometry  effects 
on  a  fully  deterministic  surface,  the  aerial  PSD  stability  analysis  evaluates 
four  contiguous  spatial  realizations  of  a  surface  that  is  nondeterministic 
but  has  underlying  spatial  structure  that  can  be  correlated. 


4.5  Spatial  analysis  results 


Introduction 


A  key  motivation  for  examining  the  spatial  attributes  of  polarimetric 
imagery  was  to  compare  the  relative  spatial  spectral  stabilities  of 
polarimetric  imagery  with  unpolarized  imagery,  recognizing  that  the  stable 
estimation  of  spatial  spectra  is  important  to  the  full  characterization  of  the 
surface  texture  and  reflection  geometry  underlying  the  spatial  radiance 
distributions  reaching  a  sensor.  Chapter  2.3  described  several  earlier 
studies  that  evaluated  the  spatial  spectra  from  unpolarized  imagery  of 
correlated  scenes.  These  results  describe  the  stablity  of  those  unpolarized 
spatial  spectral  estimates  in  comparison  with  polarimetric  estimates. 

In-field  results 


The  spatial  spectral  stabilities  of  the  three  polarimetric  components, 
I  (unpolarized  intensity),  P  (%  polarization),  and  T  (orientation  angle), 
were  analyzed  with  respect  to  a  limited  parametric  surface  exploration  of 
imaging  geometry:  the  four  scene  azimuths  and  four  solar  elevations. 

The  calibration  analysis  of  the  test  panels  in  Chapter  4.2  indicated 
that  the  lighter  of  the  two  panels  suffered  a  significant  amount  of  digital 
number  (dn)  saturation,  creating  artifacts  in  the  T  histogram  that  also 
distorted  the  calculation  of  the  T  variances;  the  saturation  would  also 
artificially  lower  P  values  since  Q  and  U  are  near  zero.  Since  the  power 
spectral  density  (PSD)  analysis  is  effectively  an  evaluation  of  the  spatial 
frequency  distribution  of  variance,  further  PSD  analysis  of  the  lighter 
panels  was  rejected  in  favor  of  the  darker  panels.  The  histograms  of  the 
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darker  panels  were  centered  within  the  linear  dn  range  and  the  resulting 
statistical  analysis  was  not  affected  by  dn  saturation. 

Figures  4.5-1  through  6  illustrate  the  sequence  of  spectral  analysis 
as  described  in  Chapter  3.4. 

The  polarimetric  component  images  in  Figures  4.5-1  through  3  (a)  are 
rotated,  fast  Fourier  transformed  (FFT),  and  the  FFT  squared  to  produce 
power  spectral  densities  (PSDs)  (in  Figures  1-3  (b))  having  the  same 
azimuth  orientation  relative  to  Scene  3. 

The  PSDs  are  combined  for  each  imaging  geometry  parameter  (scene 
azimuth  and  solar  elevation)  to  calculate  the  average  PSD  (Figures  4.5-4 
through  6  (a))  for  each  imaging  geometry  combination.  The  absolute 
differences  between  the  average  PSD  and  the  individual  PSDs  that 
formed  the  average  are  themselves  averaged  and  then  normalized  with 
respect  to  the  average  PSD  variance.  This  spatial  spectrum  is  the 
averaged  difference  PSD  (ADP).  The  advantage  of  ADP  analysis  is  that 
the  spatial  frequency  components  sum  to  a  value  that  represents  a 
percentage  of  the  total  variance  of  the  average  PSD.  This  value,  the 
integrated  difference  variance  (IDV),  is  a  single  statistic  that  indicates  the 
amount  by  which  a  group  of  PSDs  differ  from  each  other;  it  can  exceed 
100%.  The  ADP  normalization  allows  direct  comparison  of  spectra  with 
variances  expressed  in  dissimilar  units,  in  this  case,  dn  units  for  I,  % 
polarization  for  P,  and  degrees  for  T.  The  variances  of  a  group  of  images 
may  exhibit  a  large  amount  of  intra-group  variation  but  exhibit  small  IDVs 
since  it  is  the  spatial  frequency  distribution  of  relative  variation  (effectively 
the  shape  of  the  PSDs)  that  is  being  measured  and  not  absolute  variation. 
IDV,  in  sum,  is  an  unbiased  estimator  of  PSD  stability. 
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Figures  4.5-1  through  3  (a)  and  (b).  I,  P,  and  T  component  images  and  their  power 
spectral  densities  (PSDs)  for  Scene  5  &  Elevation  4.  Note  that  the  PSDs  have  been 
rotated  (in  this  case  by  +45°)  so  that  the  PSDs  for  each  scene  can  be  analyzed  with 
respect  to  a  normalized  orientation.  Also  note  that  the  underlying  spectral  components  of 
the  ‘cracked  ice’  texture  are  clearly  visible  in  all  three  PSDs  even  though  the  texture  is 
apparently  lost  in  the  P  and  T  component  images. 
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Figures  4.5-4  through  6  (a)  and  (b).  I,  P,  and  T  component  average  PSDs  and  their 
average  difference  PSDs  (ADPs)  for  Solar  Elevation  1 .  The  ADPs  are  normalized  to  the 
integrated  variance  of  the  average  PSD  and  represent  the  fractional  difference  variance 
relative  to  the  average  PSD  at  each  spatial  frequency.  The  ADP  representations  here 
were  rescaled  in  order  to  preserve  contrast  in  the  figures. 
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Because  of  several  limitations  in  the  analytic  capabilities  of  ENVI 
(discussed  in  Chapter  3.4),  fractional  IDVs  were  calculated  using  the  ADP 
histograms.  The  most  significant  limitation  was  the  quantization  (binning) 
of  the  difference  variances  during  histogram  generation  and  consequent 
aggregation  of  the  smallest  (but  most  numerous)  difference  variances  in 
the  lowest  bin.  The  information  loss  due  to  aggregation  can  be  mitigated 
by  setting  the  integration  threshhold  above  the  lowest  bin  (containing  the 
largest  expected  error)  but  the  calculated  fractional  IDV  is  still  only 
approximate  because  of  the  remaining  quantization  errors.  For  this 
reason,  several  threshholds  for  difference  variance  (dVar)  were  used  in 
the  calculation  of  fractional  IDV  in  order  to  provide  at  least  a  qualitative 
assessment  of  the  validity  of  this  approach. 

Figures  4.5-7  through  9  plot  the  fractional  IDVs  for  different  scene 
azimuths  at  dVar  threshholds  of  0.10%,  0.05%,  and  0.02%. 

Figures  4.5-10  through  12  plot  the  fractional  IDVs  for  different  solar 
elevations  at  dVar  thresholds  of  0.10%,  0.05%,  and  0.02%. 

The  IDV  values  used  in  the  plots  are  tabulated  in  Table  4-10. 

The  most  general  observation  is  that  the  fractional  IDV  values 
appear  to  preserve  relative  trends  in  the  comparison  of  the  three  dVar 
threshholds  but,  particularly,  in  the  comparison  of  the  two  lowest  dVar 
threshholds.  This  observation  provides  some  validation  for  the  fractional 
IDV  approach,  in  that  the  presence  of  increasing  quantization  errors  would 
otherwise  accumulate  during  progressive  integrations  and  disturb  the 
relative  trends. 
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With  reference  to  Figures  4.5-7  through  9,  the  IDVs  for  the  three 
polarimetric  components  I,  P,  and  T  demonstrate  a  very  clear  distinction 
of  the  relative  stabilities  of  their  PSDs  on  the  order  of  factors  2  to  3:  the 
IDVs  for  P  are  approximately  one  half  the  IDVs  of  I,  and  the  IDVs  for  T  are 
between  one  third  and  one  half  the  IDVs  of  P.  This  result  is  remarkable  in 
consideration  of  the  fact  that  P  and  T  images  are  the  product  of  ratios  that 
would,  and  do,  exhibit  increased  high  spatial  frequency  content,  similar  to 
random  noise,  and  increase  the  expected  IDV.  However,  if  this  content  is 
below  the  threshhold  for  dVar,  it  is  effectively  filtered  out  in  this  fractional 
IDV  calculation.  But,  more  robust  spectral  analysis  techniques  to  calculate 
IDV  would  otherwise  set  an  upper  spatial  frequency  threshhold  and  treat 
the  excluded  high  frequency  variance  as  noise.  The  net  result  appears  to 
be  the  same.  A  more  robust  IDV  measurement  technique  is  proposed  in 
Chapter  5. 

With  reference  to  Figure  4.5-9  alone,  the  IDVs  for  both  P  and  T  are 
within  1%  to  2%,  indicating  that  their  PSDs  are  relatively  insensitive  to 
changes  in  scene  azimuth.  While  the  2-4%  variation  in  IDVs  for  I  may  not 
be  statistically  significant,  the  apparent  trend  does  correlate  with  the 
prediction  results  of  both  Chapman  and  Irani  [1981]  and  North  [1989]: 
Scene  3  represents  the  case  where  the  dominant  spectral  components 
are  in  parallel  with  the  principal  plane,  resulting  in  maximum  errors  for 
spatial  spectrum  estimation;  Scene  5  has  the  dominant  spectral 
components  orthogonal  to  the  principal  plane,  resulting  in  minimum  errors 
for  spectrum  estimation;  and  Scenes  2  and  4  are  have  the  dominant 
components  both  45°  to  either  side  of  the  principal  plane,  resulting  in 
errors  that  are  approximately  equal  and  midway  between  the  maximum 
and  minimum  errors. 
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Figure  4.5-7  through  9,  Integrated  Difference  Variance  (IDV)  plots  by  scene  azimuth  for 
difference  variance  (dVar)  threshholds  of  0.10%,  0.05%,  and  0.02%. 


Integrated  Difference  Variance  (for  dVar  >  0.1%) 
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Integrated  Difference  Variance  (for  dVar  >  0.05%) 
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Integrated  Difference  Variance  (for  dVar  >  0.10%) 


Image  Parameter 
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Integrated  Difference  Variance  (for  dVar  >0.05%) 
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Integrated  Difference  Variance  (for  dVar  >  0.02%) 
(Combined  Scene  PSDs) 
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Figure  4.5-10  through  12.  Integrated  Difference  Variance  (IDV)  plots  by  solar  elevation 
for  difference  variance  (dVar)  threshholds  of  0.10%,  0.05%,  and  0.02%. 
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With  reference  to  Figures  4.5-10  through  12,  the  IDVs  for  the  three 
polarimetric  components  I,  P,  and  T  again  demonstrate  the  same  clear 
distinction  of  relative  stabilities  of  their  PSDs.  The  IDVs  for  both  P  and  T 
again  show  relative  insensitivity,  this  time  to  changes  in  solar  elevation. 
However,  for  increasing  solar  elevation,  the  I  parameter  IDVs  clearly 
demonstrate  a  trend  toward  increasing  instability  in  PSD  estimation.  The 
predominant  effect  in  the  imagery  is  the  presence  of  significant  shadowing 
that  increases  image  contrast  by  emphasizing  low  spatial  frequency 
changes  in  panel  elevation  and  obscuring  high  spatial  frequency  content. 
However,  the  reflected  radiance  from  the  same  shadowed  areas  is  more 
polarized  than  the  directly  illuminated  areas  since  the  source  of  incident 
radiance  comes  from  the  vicinity  of  the  antisolar  area  of  the  skydome, 
which  is  the  most  highly  polarized  part  of  the  clear  sky.  In  this  situation, 
the  higher  spatial  frequency  content  is  preserved  in  the  P  and  T  images, 
resulting  in  stabler  estimates  of  their  PSDs  relative  to  the  I  component. 

Aerial  results 


The  spatial  spectral  stabilities  of  the  three  polarimetric  components, 
I,  P,  and  T,  were  analyzed  with  respect  to  a  composite  of  four  spatially 
contiguous  realizations  (subimage  patches)  of  a  waterwave  scene 
contained  within  a  single  image. 

As  with  the  in-field  analysis  described  above.  Figures  4.5-13 
through  16  illustrate  the  sequence  of  spectral  analysis  as  described  in 
Chapter  3.4  for  the  aerial  imagery.  Figure  4.5-17  plots  the  fractional  IDVs 
at  dVar  thresholds  of  0.10%,  0.05%,  and  0.02%.  The  IDV  values  used  in 
this  plot  are  tabulated  in  Table  4-10. 
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Figures  4.5-13  through  16  (a),  (b),  and  (c).  I,  P,  T,  and  D  component  subimages  (a) 
containing  the  four  subimage  patches,  the  average  PSDs  of  the  patches  (b),  and  their 
ADPs  (c).  Again,  the  ADPs  are  normalized  to  the  integrated  variance  of  the  average  PSD 
and  represent  the  fractional  difference  variance  relative  to  the  average  PSD  at  each 
spatial  frequency.  The  ADP  representations  here  were  rescaled  in  order  to  preserve 
contrast  in  the  figures.  Also,  both  spectral  scales  were  inverted  this  time  so  that  the  larger 
magnitudes  are  represented  by  darker  pixels  and  smaller  magnitudes  by  lighter  pixels. 
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Integrated  Difference  Variance 
(Aeriai  Waterwave  Image) 


Figure  4.5-17. 

With  reference  to  Figure  4.5-17,  the  IDVs  for  the  three  polarimetric 
components  I,  P,  and  T  consistently  demonstrate  the  same  clear 
distinction  of  the  relative  stabilities  of  their  PSDs  that  was  observed  for  the 
in-field  experiment:  here  the  IDVs  for  P  are  between  one  third  and  one 
half  the  IDVs  of  I,  and  the  IDVs  forT  are  between  one  third  and  one  half 
the  IDVs  of  P. 

As  noted  in  Chapter  4.4,  the  I  image  (refer  again  to  Figure  4.5- 
13(a))  contained  significant  direct  solar  reflection  in  the  lower  portions  of 
Patches  B  and  D  that  rapidly  diminishes  in  intensity  for  pixel  values 
moving  toward  the  top  of  Patches  A  and  C;  and  the  high  contrast  would 
provide  a  more  dramatic  comparison  of  I  with  the  relative  stability  of  the  P 
and  T  spatial  spectra.  The  high  IDV  value  for  I,  in  this  case,  is  due  to  the 
four  patches  each  containing,  within  their  PSDs,  a  single  large 
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fundamental  spatial  frequency  component  that  represents  the  1 -cycle 
spatial  modulation  from  high  intensity  at  the  bottom  of  each  patch  to  low 
intensity  at  the  top.  The  differential  variance  for  the  four  PSDs  for  this  one 
spectral  component  represents  23.13%  of  the  IDV  in  the  calculation  for 
dVar  >  0.02%.  Subtraction  of  this  single  component  reduces  the  fractional 
IDV  for  I  by  almost  one  third  to  61 .22%,  providing  the  same  proportional 
scale  of  I,  P  ,  and  T  stability  for  the  aerial  image  as  is  observed  for  the  in¬ 
field  images. 

The  ADP  was  also  calculated  for  the  difference  (D)  subimage 
patches  to  give  an  indication  of  the  spatial  spectral  stability  of  the 
systematic  error  in  the  uncalibrated  image.  While  D  is  highly  correlated 
with  the  other  three  spectra,  it  also  has  the  highest  relative  stability  of  the 
four  ADPs,  both  factors  suggesting  that  the  increased  aerial  IDV  values 
for  I,  P,  and  T  relative  to  the  corresponding  in-field  IDV  values  are  due  to 
the  calibration  errors  associated  with  D.  The  relative  stability  of  D  also 
suggests  that  the  increase  in  the  aerial  IDV  values  created  by  the  non¬ 
zero  presence  of  D  is  in  relatively  constant  proportion  to  the  calibrated  in¬ 
field  IDV  values.  This  appears  to  be  the  case:  the  uncalibrated  aerial  IDV 
values  are  larger  than  their  corresponding  calibrated  in-field  IDV  values  by 
a  constant  factor  of  about  two.  The  inverse  of  this  statement  is  that,  based 
on  the  limited  results  here,  sensitometric  calibration  of  this  image  data 
appears  to  roughly  double  the  stability  of  the  spectral  estimates  compared 
with  uncalibrated  image  data.  Enhanced  quantitative  analysis  and 
validation  of  these  limited  observations  is  recommended  in  Chapter  5. 
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Table  4-10 
IDV  Summary  Report 


In-field  results  by  scene  (for  dVar  >  0.10%) 


Parameter  Scene  2  Scene  3  Scene  4  Scene  5 


1 

11.81 

14.32 

11.50 

7.76 

p 

1.35 

1.58 

1.65 

1.80 

T 

0.00 

0.00 

0.56 

0.25 

(for  dVar  >  0.05%) 

1 

18.94 

23.42 

17.68 

16.46 

P 

5.50 

4.47 

6.18 

4.71 

T 

0.73 

1.01 

1.17 

0.92 

(for  dVar  >  0.02%) 

1 

28.85 

32.56 

26.51 

24.17 

P 

16.10 

12.96 

14.45 

12.77 

T 

5.27 

6.33 

5.33 

5.01 

In-field  results  by  solar  elevation  (for  dVar 

>0.10%) 

Parameter  Sol  Elev  1 

Sol  Elev  2 

Sol  Elev  3 

Sol  Elev  4 

1 

10.70 

12.78 

14.04 

25.36 

P 

0.69 

3.25 

2.51 

1.83 

T 

0.30 

0.00 

0.00 

0.44 

(for  dVar  >  0.05%) 

1 

20.05 

22.97 

25.86 

39.65 

P 

6.76 

8.03 

8.95 

6.06 

T 

1.00 

1.28 

0.92 

1.24 

(for  dVar  >  0.02%) 

1 

30.24 

33.06 

35.93 

54.21 

P 

18.24 

16.53 

19.70 

17.03 

T 

4.72 

7.59 

6.59 

5.90 

Aerial  results  by  dVar 

Parameter  >0.10% 

>  0.05% 

>  0.02% 

1 

82.18 

84.35 

84.35 

P 

15.20 

23.55 

33.01 

T 

5.11 

9.08 

16.29 

D 

3.76 

5.75 

13.14 
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5.0  Conclusions  and  recommendations 


Conclusions 


The  results  of  this  study  have  produced  five  main  conclusions: 

1)  Sensitometric  co-calibration  of  the  characteristic  curves  for  the 
input  images  (i.e.,  the  photo  quads)  is  important  for  the  accurate 
calculation  of  the  Stokes  parameters  and  derivatives.  In  particular,  the 
non-linear  transformation  of  intensity  is  shown  to  distort  the  calculation  of 
T  and  its  variance  for  the  extreme  case  of  saturation.  Saturation  at  the 
extrema  should  be  avoided;  linear  transformation  over  a  wide  range  of 
intensity  should  be  sought. 

2)  The  polarimetric  difference  image  (D)  provides  a  direct 
measurement  of  sensitometric  co-calibration  error  and  a  potential  tool  for 
mitigating  systematic  error. 

3)  For  the  reflecting  surfaces  that  were  imaged  under  clear  sky 
conditions,  the  histogram  distributions  of  T  provide  an  indirect  measure  of 
the  azimuthal  difference  between  the  polarizer  reference  plane  and  the 
principal  plane  containing  the  sun.  This  indirect  measurement  provides 
the  potential  capability  to  calculate  the  the  azimuth  of  either  plane  when 
only  one  is  known,  or  to  calibrate  the  system  when  the  azimuths  of  both 
are  known.  The  measurements  acquired  within  the  current  study  indicate 
that  the  calibrated  polarimeter  values  of  T  were  in  good  agreement  with 
independent  calculations. 
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4)  The  calibrated  polarimetric  components  P  and  T  provide  a  more 
stable  estimate  of  spatial  PSD  distribution  in  relation  to  I.  Within  the 
limited  exploration  of  parameters,  P  demonstrates  PSDs  that  have  twice 
the  stability  of  1;  and  T  demonstrates  PSDs  that  have  two  to  three  times 
the  stability  of  P. 

5)  Within  the  limited  exploration  of  parameters,  the  stability  of  the 
PSDs  of  the  calibrated  polarimetric  components  P  and  T  demonstrate  a 
relative  insensitivity  to  imaging  geometry  compared  with  I. 

In  sum,  the  results  of  this  study  indicate  that  the  calibrated 
polarimetric  image  set  (I,  P,  T)  has  the  potential  to  provide  a  significant 
advantage  over  the  unpolarized  image  component  (I)  alone  in  the 
detection,  discrimination,  and  classification  of  spatial  content  through  the 
relative  stability  and  geometric  invariance  in  the  estimation  of  spatial 
power  spectral  density. 


Recommendations 

The  list  of  recommendations  follows  the  sequential  analysis  of  the 
full  polarimetric  imaging  chain,  beginning  with  the  illumination  source  and 
ending  with  interpretation  of  the  final  image  products. 

Illumination  source  /  radiation  path  analysis 

1)  Further  imaging  studies  should  provide  detailed  insight  into  the 
full  range  of  natural  polarized  radiance  conditions  that  can  illuminate  both 
natural  and  artificial  terrestrial  scene  contents.  In  particular,  optimizing 
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illumination  conditions  for  polarimetric  remote  sensing  should  be  explored. 
A  full  parametric  surface  exploration  of  imaging  geometry  combinations 
under  various  natural  polarized  radiance  conditions  is  recommended. 

Polarized  target  characterization 

2)  Enhanced  imaging  studies  should  emphasize  the  character¬ 
ization  of  polarized  bidirectional  reflectance  distributions  (BRDF)  for 
terrestrial  features  under  natural  polarized  radiance  conditions.  In 
particular,  the  polarized  spatial  content  analysis  of  large  spatially 
correlated  and  uncorrelated  natural  features  should  be  investigated  for 
optimizing  imaging  geometries. 

3)  The  eight  targets  of  opportunity  that  were  imaged  during  the 
aerial  collection  exist  as  a  point  of  departure  for  an  initial  color  remote 
sensing  assessment. 

Polarimetric  imaging  sensor  enhancements 

4)  Other  linear  polarizers,  for  example  Polaroid  HN38S,  should  be 
evaluated  for  potentially  improved  remote  sensing  performance  in  the 
visible  spectrum;  an  investigation  of  infrared  and  multispectral  polarimetric 
imaging  strategies  should  also  be  considered,  for  instance,  in  a  ganged, 
multi-camera,  multi-lens  camera  system  that  could  emulate  the  POLDER 
instrument  but  with  simultaneous,  rather  than  sequential,  image  collection. 

5)  If  a  different  camera  can  be  employed,  high  quality  lenses 
should  be  incorporated  in  a  multi-lens,  mapping  camera  type  system.  The 
system  should,  as  a  minimum,  include  positive  shutter  control  and 
motorized  film  advance  (features  that  are  lacking  in  the  current  system). 
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6)  Combinations  off-stop,  exposure  time,  film  speed,  and 
polarizing  filter  density  should  be  evaluated  for  their  ability  to  provide  the 
largest  linear  range  of  exposure  (i.e.,  exposure  latitude)  for  polarimetric 
remote  sensing  in  support  of  the  calibrated  transformation  of  exposure  to 
intensity. 

7)  If  possible,  detector  arrays  should  be  considered  as  a 
replacement  for  film  at  the  four  focal  planes.  The  most  significant  benefit 
of  current  solid  state  charge-coupled  device  (CCD)  arrays  is  their 
excellent  linearity  in  output  over  a  very  large  range  of  input  intensity.  This 
benefit  obviates  the  requirement  for  sensitometric  procedures  except  for 
initial  system  calibration  and  absolute  calibration  maintainance.  The  many 
other  benefits  include  direct  digitization,  contrast  control,  direct  radiometric 
co-calibration  control,  and  direct  co-registration  control.  The  ultimate 
advantage  of  this  replacement  is  the  potential  for  near-real-time  creation 
of  the  polarimetric  image  set  since  the  arrays  could  port  their  image  data 
directly  into  a  digital  computer  and  Stokes  parameter  processing  could 
automatically  initiate  upon  receipt  of  quad  image  data. 

Polarimetric  image  processing  enhancements 

8)  Additional  film  digitization  techniques  should  be  evaluated  for 
their  ability  to  provide  faithful  reproduction  of  digital  number  values  that 
linearly  and  directly  correlate  with  relative  intensity  values.  For  instance,  a 
dedicated  film  scanner  would  provide  the  researcher  with  full  control  over 
the  digitization  process,  in  particular,  control  of  contrast  and  color  balance. 
Several  35-mm  slide  scanners  exist  in  the  market  that  have  excellent 
illumination  and  thermal  stability  because  they  use  light  emitting  diode 
(LED)  technology;  one  unit  also  has  high  spatial  resolution  (9.4-micron 
spot  size)  and  fits  within  a  standard  5.25-inch  computer  drive  bay.  The 
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flexibility  of  these  scanning  systems  also  allows  consideration  of  color 
reversal  film  (e.g.,  slide  transparencies) 

9)  One  significant  improvement  for  future  Stokes  vector  processing 
would  be  the  incorporation  of  a  fully  automated  or  semi-automated  co¬ 
registration  capability  using  cross-correlation.  This  enhancement  would 
remove  the  operator  as  a  source  of  registration  error  through  elimination 
of  the  operator’s  tedious  and  subjective  selection  of  registration  points  that 
are  common  within  the  four  image  quads.  For  example,  co-registration 
within  the  current  study  involved  the  selection  of  eight  common  points  for 
each  of  the  four  quads,  or  32  total  points  for  each  Stokes  image  set. 


Figure  5.0-1  illustrates 
an  example  of  cross¬ 
correlation.  Image  (b) 
is  translated  relative  to 
image  (a),  the  product 
of  their  pixel  values 
are  summed  for  each 
translation,  and  the 
sum  reported  as  the 
pixel  value  of  image 
(c). 


Figure  5.0-1  (a),  (b),  and  (c).  Example  of  cross-correlation 
[from  Gonzalez  &  Wintz,  1987]. 


The  enhancement  would  involve  the  cross-correlation  of  small  (e.g.,  20  x 
20  pixel)  unregistered  subimages,  both  containing  the  same  approximate 
in-scene  features,  to  detect  the  common  centroid  coordinates  as  the 
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highest  peak  in  the  cross-correlation  through  the  calculation  of  the 
normalized  cross-correlation  coefficient: 

r-(vv\-  +  +  I  Equation  5-1 

^ .dC^  ^  Jl  y  j - -  - -  ^ 

where  f(x,y)  and  g(x,y)  are  the  two  subimages.  C(X,Y)  attains  a  maximum 
value  when  the  two  subimages  match  each  other,  as  illustrated  in  Figure 
5.0-1  (c). 

The  main  simplifying  assumption  is  that  the  co-registration  of  the  image 
quads  is  predominantly  a  translation  of  coordinates  rather  than  a  rotation 
or  scale  change,  since  the  four  optic  axes  and  four  focal  plane  axes  are  all 
fixed  approximately  in  parallel  with  each  other;  therefore,  the  cross¬ 
correlation  effectively  evaluates  all  translation  combinations  of  subimage 
coordinates  until  the  spatial  distributions  of  the  two  subimages  achieve  a 
match.  Multiple  centroid  coordinates  can  be  automatically  generated  until 
the  number  is  sufficient  for  the  polynominal  fit  of  the  desired  registration. 

An  evaluation  of  this  co-registration  approach  is  recommended. 

10)  A  second  significant  improvement  for  future  Stokes  vector 
processing  would  be  the  incorporation  of  a  fully  automated  or  semi- 
automated  sensitometric  cross-calibration  capability.  This  capability  would 
take  advantage  of  the  ready  availability  of  the  difference  (D)  image  as  a 
derivative  analytic  product  for  co-calibration  refinement.  The  key 
simplifying  assumption  is  that  the  recorded  exposures  of  the  four  image 
quads  correlate  with  intensity  through  a  single  characteristic  curve  (an 
assumption  that  should  have  held  but  did  not  hold  for  the  current  study). 
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The  residual  error  in  the  D  image  should  then,  for  the  most  part,  represent 
the  effect  of  non-linearity  in  the  characteristic  curve  at  the  extrema  and,  to 
a  lesser  degree,  represent  the  effect  of  other  system  errors  such  as  pixel 
mis-registrations  from  parallax,  lens  and  filter  vignetting,  and  etc. 

This  enhancement  would  essentially  involve  a  pixel-by-pixel  analysis  of 
the  four  quad  images  and  the  difference  (D)  image  with  the  intent  of 
minimizing  D  by  equalizing  the  sum  of  the  two  polarized  pairs,  i.e.,  making 
r  equal  to  I”.  The  simplest  case  for  the  analysis  logic  would  first  evaluate 
the  two  polarized  pairs  to  determine  which  of  the  quads  have  dn  values 
that  are  nearest  the  extrema  of  the  characteristic  curve  and  suffer  most 
from  nonlinearity  in  transformation.  If  three  of  the  four  dn  values  fall  within 
the  linear  region,  the  value  of  D  would  be  applied  to  the  fourth  value  so 
that  r  =  I”  and  D  =  0.  A  more  sophisticated  approach  would  be  to  estimate 
the  amount  of  nonlinearity  that  each  pixel  value  in  each  quad  may  suffer, 
calculate  the  histogram  distribution  of  these  estimates,  and  then  use  an 
iterative  least  squares  optimization  approach  to  distribute  a  weighted 
fraction  of  D  to  each  of  the  quads  until  the  difference  and  the  histogram 
both  approach  acceptable  minima. 

An  evaluation  of  this  calibration  approach  is  recommended. 

Polarimetric  image  analysis  enhancements 

11)  Enhanced  imaging  studies  should  consider  more  robust  PSD 
stability  measurements,  using  a  full  parametric  surface  exploration  of 
imaging  geometries  for  terrestrial  features  under  natural  polarized 
radiance  conditions.  In  particular,  PSD  stability  analysis  of  large  spatially 
correlated  and  uncorrelated  natural  features  should  be  investigated  for 


150 


optimizing  imaging  geometries.  This  effort  should  hopefully  provide  an 
independent  validation  of  the  limited  IDV  results  in  this  study. 

Enhanced  PSD  and  PSD  stability  measurements  would  include  the 
capability  to  execute  the  double  integration  of  the  two-dimensional  ADP 
surface,  instead  of  the  approximate  histogram  method  employed  in  this 
study.  The  most  significant  benefit  of  this  capability  is  the  option  to 
selectively  calculate  IDV  as  a  function  of  one  variable,  spatial  frequency  or 
azimuth  (through  a  partial  integration  over  the  other  variable),  and  then 
evaluate  the  stability  of  the  PSD  with  respect  to  these  spatial  parameters: 


0=0 

Ir 

'^max 

s(e)=i;s.(e) 

A=0 


Equation  5-2 


Equation  5-3 


Figure  5.0-2  illustrates  examples  of  the  one-dimensional  spatial  spectra  in 
comparison  with  the  original  image  and  its  two-dimensional  spectrum. 

An  evaluation  of  this  spatial  analysis  enhancement  is  recommended. 

12)  Finally,  a  calibrated  polarimetric  imaging  model  should  be 
validated  with  respect  to  the  performance  of  a  calibrated  imaging 
polarimeter.  The  results  of  Chapter  4.3  provide  some  qualitative  insight 
into  the  spatial  performance  of  a  limited  model  imaging  a  specific  feature 
within  a  limited  range  of  imaging  geometries  under  clear  sky  conditions; 
but  a  robust  model  needs  to  provide  the  quantitative  correlation  that  can 
predict  system  performance  for  the  more  general  case  of  polarimetric 
remote  sensing  of  various  terrestrial  features  collected  over  a  broad  range 
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Figure  5.0-2  (a),  (b),  (c),  and  (d).  Image  (a),  its  two-dimensional  spatial  spectrum 
S(spatial  frequency,  azimuth)  (b),  its  one-dimensional  spectrum  S(spatial  frequency)  (c), 
and  its  one-dimensional  spectrum  S(azimuth)  (d)  [from  Gonzalez  &  Wintz,  1987]. 


of  imaging  geometries  and  illumination  conditions,  many  of  which  may  not 
be  ideal. 
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Appendix  A  Source  code  for  the  synthetic  imaging  model 


PROGRAM  SKYIMAGE 

* 

*  Creates  synthetic  polarized  reflected  skydome  intensity  images  as  functions 

*  of  variables  Beta  &  Alpha  (dome  mirror  slope  angles)  &  fixed  constants  Theta 

*  (sensor  declination),  Phi  (sensor  azimuth:  fixed  at  180  deg),  ThetO  (sun 

*  declination)  &  PhiO  (sun  azimuth).  Other  program  constants: 

* 

*  FD  (focal  distance  from  mirror) 

*  MRD  (radius  of  dome  mirror) 

*  DSD  (pixel  detector  sampling  distance) 

*  Nreal  &  Nimag  (complex  surface  reflection  coefficients  of  mirror) 

*  K_1  &  KJZ  (transmission  coefficients  for  non-ideal  linear  polarizer) 

*  Depol  (depolarization  coefficient  for  non-ideal  linear  polarizer) 

*  Gamma  &  Bias  (effective  sensor  exposure  gamma  &  offset) 

*  DN_max  (maximum  digital  number  for  pixel  intensity) 

* 

*  [Intensity  arrays  span  Beta  angle  values  <=  45  degrees,  to  correspond  with 

*  the  upper  bound  for  surface  slope  that  reflects  modeled  skydome  radiance] 

* 

*  Creates  derivative  Stokes  parameter  images  from  the  synthetic  polarized 

*  reflected  skydome  intensity  images. 

* 

IMPLICIT  NONE 

INTEGER  l,J.N,P,Q 

PARAMETER  (N=64,P=-N/2,Q=N/2-1) 

REAL  ISRF  [FAR](P:Q,P:Q),QSRF  [FAR](P:Q,P;Q),USRF  [FAR](P:Q.P:Q) 
REAL  PSRF  [FAR](P:Q,P:Q),TSRF  [FAR](P:Q,P:Q) 

REAL  1000  [FAR](P;Q,P:Q),I045  [FAR](P:Q,P:Q).1090  [FAR](P:Q,P;Q) 
REAL  1135  [FAR](P:Q,P;Q) 

REAL  DSD, FD,FL, MRD, PHIO, PI, THETA.THET0, NREAL, NIMAG, N_MAG 
REALPARM(1:4,1:4),K_1,K_2,POL_MAX,DEPOL,GAMMA,BIAS,DN_MAX 
CHARACTER*24  FNAME,FNAME1  ,FNAME2,FNAME3,FNAME4,FNAME5 
CHARACTER*24FNAME6,FNAME7,FNAME8,FNAME9 

* 

*  Define  constants 

* 

PI  =  3.141593 

* 

*  Input  parameters  from  external  file 

* 

WR1TE(*,*)  ’  Enter  filename  of  input  parameters  ' 

READ(*,300)  FNAME 

WRITE  (*,*)  ’  Reading  parameter  file . ' 

OPEN(1  ,FILE=FNAME,ERR=5) 

DO  l=1,4 

READ(1,*)  (PARM(I,J),J=1,4) 

ENDDO 

CLOSE(I) 

* 

WRITE(*,*)  ’  ’ 

DO  1=1,4 

WRITE(*.*)  (PARM(I,J), J=1 ,4) 

ENDDO 

PAUSE 

* 

*  Assign  variable  names  to  input  parmaters 

* 

FD  =  PARM(I.I) 

FL  =  PARM(1,2) 

MRD  =  PARM(1,3) 

DSD  =  PARM(1,4) 

POL_MAX  =  PARM(2.1) 

DEPOL  =  PARM(2,2) 
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GAMMA  =  PARM(2,3) 

BIAS  =  PARM(2.4) 

THETA  =  PARM(3,1) 

PHIO  =  PARM(3,2) 

THETO  =  PARM(3,3) 

DN_MAX  =  PARM(3,4) 

NREAL=PARM(4,1) 

NIMAG  =  PARM(4,2) 

KJ  =  PARM(4.3) 

K_2  =  PARM(4,4) 

* 

THETA  =  THETA*(PI/180.0) 

PHIO  =  PHI0*(P1/1 80.0) 

THETO  =  THET0*(PI/180.0) 

N„MAG  =  SQRT(NREAL**2.0  +  NIMAG**2.0) 

it 

*  Call  subroutine  SKYSURF 

ic 

CALL  SKYSURF  (PHIO, THETO, THETA, N_MAG.FD,FL,MRD, DSD, POL_MAX,DEPOL, 
+  GAMMA,BIAS,DN_MAX,K_1  ,K_2,P,Q,ISRF,QSRF,USRF, 

+  PSRF,TSRF,IOOO,1045,I090,I135) 

ic 

WRITER, *)  • ' 

WRITE(*,*) '  SKYSURF  completed  successfully' 

* 

WRITE(*  *) ' ' 

WRITE(*  *) '  Enter  filename  w/  drive  letter  &  directory  for 
+  '  ISRF  data' 

*  READ(*,300)  FNAME1 

FNAME1  =  'ISRF.DAT' 

WRITER, *)  ’  ’ 

WRITE(*,*)  ’  Enter  filename  w/  drive  letter  &  directory  for 
+  '  QSRF  data' 

*  READ(*,300)  FNAME2 

FNAME2  =  'QSRF.  DAT' 

WRITEf,*) ' ' 

WRITE(*,*) '  Enter  filename  w/  drive  letter  &  directory  for 
+  ’  USRF  data' 

*  READ(*,300)  FNAME3 

FNAME3  =  'USRF.DAT' 

WRITEf,*) ' ' 

WRITE(*,*) '  Enter  filename  w/  drive  letter  &  directory  for 
+  ’  PSRF  data’ 

*  READ(*.300)  FNAME4 

FNAME4  =  'PSRF.DAT' 

WRITER  *) ' ' 

WRITE(*  *) '  Enter  filename  w/  drive  letter  &  directory  for 
+  ’  TSRF  data’ 

*  READr,300)  FNAME5 

FNAME5  =  ’TSRF.DAT' 

WRITE(*,*)  ’ ' 

WRITER,*) '  Enter  filename  w/  drive  letter  &  directory  for 
+  ’  lOOO  data' 

*  READ(*,300)  FNAME6 

FNAME6  =  ’I000.DAT’ 

WRITE(*  *) '  ’ 

WRITE(*,*)  ’  Enter  filename  w/  drive  letter  &  directory  for 
+  ’  I045  data’ 

*  READ(*,300)  FNAME7 

FNAME7  =  ’1045.DAT 
WRITER, *) ' ' 

WRITE(*,*) '  Enter  filename  w/  drive  letter  &  directory  for 
+  ’  I090  data' 

*  READ(*,300)  FNAME8 

FNAME8  =  ’lOQO.DAT’ 

WRITER, *) '  ’ 

WRITER,*)  ’  Enter  filename  w/  drive  letter  &  directory  for 
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+  ’1135  data' 

READ(*.300)  FNAME9 

FNAME9  = ’I135.DAT’ 

WR1TE(*,*)  ’  ’ 

OPEN(2,FILE=FNAME1,FORM=’FORMATTED’.ERR=5) 
DO  l=P,Q 

WRITE(2,*)  (ISRF(1,J),J=P.Q) 

ENDDO 

CLOSE(2) 

WRITEf,*) '  '.FNAMEi;  has  been  output' 

OPEN(2,FILE=FNAME2,FORM='FORMATTED',ERR=5) 
DO  l=P,Q 

WRITE(2.*)  (QSRF(I.J),J=P,Q) 

ENDDO 

CLOSE(2) 

WRITER,*)  ’  ',FNAME2;  has  been  output’ 

OPEN(2,FILE=FNAME3,FORM='FORMATTED',ERR=5) 
DO  l=P.Q 

WRITE(2,*)  (USRF(I.J),J=P,Q) 

ENDDO 

CLOSE(2) 

WRITEf,*) '  ',FNAME3,’  has  been  output’ 

OPEN(2,FILE=FNAME4.FORM='FORMATTED',ERR=5) 
DO  l=P,Q 

WRITE(2  *)  (PSRF(I,J),J=P.Q) 

ENDDO 

CLOSE(2) 

WRITE(*,*)  ’  ’,FNAME4,’  has  been  output’ 

OPEN(2,FILE=FNAME5,FORM='FORMATTED’,ERR=5) 
DO  l=P,Q 

WRITE(2,*)  (TSRF(U),J=P,Q) 

ENDDO 

CLOSE(2) 

WRITE(*,*) '  ',FNAME5;  has  been  output' 

OPEN(2,FlLE=FNAME6,FORM=’FORMATTED’,ERR=5) 
DO  1=P,Q 

WRITE(2 ,*)  (I000(I,J),J=P,Q) 

ENDDO 

CLOSE(2) 

WRITEf  *) '  '.FNAMES;  has  been  output’ 

0PEN(2.F1LE=FNAME7,F0RM='F0RMATTED',ERR=5) 
DO  l=P,Q 

WRITE(2,*)  (1045(1, J).J=P,Q) 

ENDDO 

CLOSE(2) 

WR!TE(*,*)  ’  ',FNAME7,'  has  been  output' 

OPEN(2,FILE=FNAME8,FORI\/l=’FORMATTED',ERR=5) 
DO  l=P,Q 

WR1TE(2,*)  (I090(I.J).J=P,Q) 

ENDDO 

CLOSE(2) 

WRITE(*  *)  ’  ’,FNAME8,’  has  been  output’ 

OPEN(2,FILE=FNAME9,FORM=’FORMATTED',ERR=5) 
DO  l=P.Q 

WRITE(2,*)  (1135(1, J),J=P,Q) 

ENDDO 

CLOSE(2) 

WRITE(*.*) '  ’,FNAME9,'  has  been  output' 


G0T0  1 

5  WRITEf,*)  'Bad  file  opening’ 

300  FORMAT  (A) 

1  END 

* 

SUBROUTINE  SKYSURF  (PHIO.THETO, THETA, N_MAG,FD,FL,MRD, DSD, POL_MAX, 
+  DEPOL,GAMMA,BIAS,DN_MAX,K_1  ,K_2,P,Q, 

+  ISRF,QSRF,USRF,PSRF,TSRF, 1000,1045,1090,11 35) 

* 

IMPLICIT  NONE 
INTEGER  l,J,P,Q 

REAL  ALPH,BETA.COSOMEGA,DEGC,THET,THETA,HALFPI,PI,QTRPI 
REAL  ITMP,IX,IY,IZ,MAGI,NX,NY,NZ,OMEGA,QTMP,RX,RY,RZ,UTMP 
REAL  ISRF  [FAR](P;Q,P:Q),QSRF  [FAR](P:Q,P:Q),USRF  [FAR](P:Q,P:Q) 

REAL  PSRF  [FAR](P:Q,P:Q),TSRF  [FAR](P:Q,P:Q) 

REAL  lOOO  [FAR](P:Q,P:Q),I045  [FAR](P:Q,P:Q),I090  [FAR](P:Q,P:Q) 

REAL  II 35  [FAR](P:Q,P:Q) 

REAL  A,B,BIAS,C,COSMU,COSMU2,COSNU,DEPOL,GAMMA,L,LREF,M,MU,NU,PHI 
REAL  PHI0,POL_MAX,PSI,PX,PY,PZ,QX,QY,QZ,SX,SY,SZ,THET0 
REALATT,MDIAG,MDIAQ,MOFFD,N_MAG,RHO.RPAR,RPRP,TTMP 
REAL  DN_MAX,DSD,FD,FL,IMAX1,IMAX2,K_1,K_2,IRAD,MHT,MRD,MRAD 

* 

*  Define  constants 

* 

PI  =  3.141593 
DEGC=  PI/180.0 
HALFPI  =  PI/2.0 
QTRPl  =  PI/4.0 

*  Initiate  variables 

* 

IMAX1  =  0.0 
IMAX2  =  0.0 


*  Calculate  LREF  (zenith  radiance  value)  relative  to  direct  sun  radiance 

*  (set  at  100.0)  based  on  CIE  sky  radiance  model  (after  Hopkinson) 

* 

COSMU=  1.0 
MU  =  0.0 

A  =  0.91+10.0*EXP(-3.0*MU)+0.45*(COSMU**2.0) 

* 

IF  (THETO.LT.HALFPI)  THEN 
B  =  1.0-EXP(»0.32/COS(THET0)) 

ELSEIF  (THETO.EQ.HALFPI)  THEN 
B  =  1.0 
ELSE 

WRITE(*,*) ' ' 

WRITE(*  *)  'Solar  declination  indicates  sun  is  below  horizon' 
ENDIF 

* 

C  =  0.274*(0.91+10.0*EXP(-3.0*THET0)+0.45*(COS(THET0)**2.0)) 
LREF  =  (100.0*C)/(A*B) 

* 

*  Analyze  reflected  radiance  unit  vector  R  into  xyz  components 

RX  =  SIN(THETA)*COS(PI) 

RY  =  SIN(THETA)*SIN(PI) 

RZ  =  COS(THETA) 

* 

*  Initiate  double  loop  to  generate  2-D  image  arrays 

* 

DO  I  =  P+1, Q 
DO  J  =  P+1.Q 

* 

*  Calculate  pixel  coordinates 
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*  [Row  1/L  specifies  Y-magnitude;  Column  J/M  specifies  X-magnitude] 

★ 

L  =  FLOAT(l)*DSD 
M  =  FLOAT(J)*DSD 

* 

*  Calculate  slope  angle  (BETA)  and  azimuth  angle  (ALPH)  of  mirror  facet 

*  imaged  by  pixel  (l,J) 

* 

*  [Calculate  radius  of  mirror  image  in  image  coordinates] 

*  [Calculate  radius  of  mirror  in  cartesian  (xy)  coordinates:  x/X  =  fl/FD] 

*  [Calculate  height  of  mirror  facet  in  cartesian  (z)  coordinates] 

* 

IRAD  =  SQRT((L**2.0)+(M**2.0)) 

MRAD  =  (FD*IRAD)/FL 

MHT  =  SQRT(IVIRD**2.0  -  MRAD**2.0) 


*  [Calculate  BETA] 

* 

IF  (MHT.GT.0.0)  THEN 

BETA  =  ATAN(MRAD/MHT) 

ELSE 

BETA  =  0.0 

ENDIF 

* 

*  Set  image  arrays  to  zero  if  incident  vector  arrives  from  below  the  horizon 

*  [i.e.,  if  2*OMEGA  >  90  deg;  or  if  THET  >  90  deg] 

* 

IF  (BETA.GT.QTRPI)  THEN 
ISRF(l,J)  =  0.0 
QSRF(I,J)  =  0.0 
USRF(I,J)  =  0.0 
PSRF(I.J)  =  0.0 
TSRF(l,J)  =  0.0 
I000(I,J)  =  0.0 
1045(1, J)  =  0.0 
1090(1, J)  =  0.0 
I135(I,J)  =  0.0 
GOTO  4 

ENDIF 

■k 

*  [Calculate  ALPH] 

* 

IF  (BETA.GT.0.0)  THEN 

ALPH  =  ATAN2(L,M) 

ELSE 

ALPH  =  0.0 

ENDIF 

*  Analyze  surface  normal  unit  vector  N  into  xyz  components 

* 

NX  =  SIN(BETA)*COS(ALPH) 

NY  =  SIN(BETA)*SIN(ALPH) 

NZ  =  COS(BETA) 

*  Dot  product  of  R  &  N  =  |R||N|cos(omega)  =  (1 .0)  cos(omega) 

* 

COSOMEGA  =  RX*NX  +  RY*NY  +  RZ*NZ 
OMEGA  =  ACOS(COSOMEGA) 


*  Specify  xyz  components  of  I,  the  incident  skydome  radiance  vector 

IX  =  (NX*2.0*COSOMEGA)  -  RX 
lY  =  (NY*2.0*COSOMEGA)  -  RY 
IZ  =  (NZ*2.0*COSOMEGA)  -  RZ 

* 

*  Normalize  I  ->  |l|  =  1.0 
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MAGI  =  SQRT(IX**2.0  +  1Y**2.0  +  IZ**2.0) 

IX  =  IX/MAGI 
lY  =  lY/MAGI 
IZ  =  IZ/MAGI 

*  Calculate  incident  skydome  angles  THEY  &  PHI 

* 

THEY  =  ACOS(IZ) 

ie 

*  Set  image  arrays  to  zero  if  incident  vector  arrives  from  below  the  horizon 

*  [i.e.,  if  2*OMEGA  >  90  deg;  or  If  YHEY  >  90  deg] 

* 

IF  (YHEY.GY.HALFP!)  YHEN 
ISRF(I,J)  =  0.0 
QSRF(I,J)  =  0.0 
USRF(I,J)  =  0.0 
PSRF(I,J)  =  0.0 
YSRF(U)  =  0.0 
1000(1, J)  =  0.0 
1045(1, J)  =  0.0 
1090(1, J)  =  0.0 
I135(I,J)  =  0.0 
G0Y0  4 

ENDIF 

* 

IF  (YHEY.GY.0.0)  YHEN 

PHI=AYAN2(IY,IX) 

ELSE 

PHI  =  0.0 

ENDIF 

* 

*  Calculate  radiance  magnitude  of  incident  skydome  vector  I 

*  [CIE  sky  radiance  model  (after  Hopkinson)] 

* 

1  COSMU  =  COS(YHEYO)*COS(YHEY)+SIN(YHEYO)*SIN(YHEY)*COS(PHIO-PHI) 
MU  =  ACOS(COSMU) 

A  =  0.91+10.0*EXP(-3.0*MU)+0.45*(COSMU**2.0) 

★ 

IF  (YHEY.LY.HALFPI)  YHEN 

B  =  1.0-EXP(-0.32/COS(YHEY)) 

ELSE 

B=  1.0 

ENDIF 

2  IF(MU.EQ,0.0)YHEN 

PSI  =  0.0 
COSNU  =  1.0 
NU  =  0.0 

ELSE 

COSMU2  =  COSMU**2.0 

PSI  =  POL_MAX*((1.0-COSMU2)/(1.0+COSMU2)) 

*  P  =  Polarization  Vector  (vector  orientation  of  the  polarization  ellipse 

*  =  cross  product  of  Sun  Vector  &  Sky  Element  Vector  =  S  x  Q 

* 

*  Analyze  unit  vector  Q  &  S  into  xyz  components 

*[Q  =  |I|] 

* 

QX  =  S!N(YHEY)*COS(PHI) 

QY  =  SIN(YHEY)*SIN(PHI) 

QZ  =  COS(YHEY) 

* 

=  SIN(YHEY0)*COS(PHI0) 

SY  =  SIN(YHEY0)*SIN(PHI0) 

SZ  =  COS(YHEYO) 

★ 

*  Calculate  the  cross  product,  P  =  S  x  Q 
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PX  =  SY*QZ  -  SZ*QY 
PY  =  SZ*QX  -  SX*QZ 
PZ  =  SX*QY  -  SY*QX 

* 

*  Calculate  orientation  angle  NU  of  the  polarization  ellipse 

*  Rem:  Y-axis  runs  +  to  -  for  quadrant  1  &  II  calculations 

* 

IF  (ABS(PY).GT.O.O)  THEN 
NU  =  ATAN(PX/PY) 

ELSE 

NU  =  0.0 

ENDIF 

* 

COSNU  =  COS(NU) 

ENDIF 

* 

1SRF(I,J)  =  LREF*A*B/C 

QSRF(U)  =  1SRF(I,J)*PS1*COS{2.0*NU) 

USRF(I,J)  =  ISRF(l,JrPSI*SIN(2.0*NU) 

PSRF(I,J)  =  PSI 
TSRF(I,J)  =  NU 

* 

*  Rotate  Stokes  vector  to  reference  plane  of  reflecting  facet 

QTMP  =  COS(-ALPHrQSRF(l,J)+SIN(-ALPH)*USRF(l,J) 

UTMP  =  -SIN(-ALPHrQSRF(I.J)+COS(-ALPH)*USRF(l,J) 
QSRF(I,J)  =  QTMP 
USRF(I.J)  =  UTMP 

* 

*  Calculate  perpendicular  &  parallel  reflection  coefficients  for 

*  incident  angle  OMEGA,  refraction  angle  RHO,  &  refractive  index  N_MAG 

* 

IF  (ABS(OMEGA).EQ.O.O)  OMEGA  =  0.00000001 

* 

RHO  =  ASIN(SIN(OMEGA)/N_MAG) 

rPAR  =  -TAN(OMEGA-RHO)/TAN(OMEGA+RHO) 

RPRP  =  -SIN(OMEGA-RHO)/SIN(OMEGA+RHO) 


*  Specify  Mueller  matrix  for  reflecting  dielectric  surface 

*  [MDIAG  =  Mil  =  M22,  MOFFD  =  M12  =  M21,  MDIAQ  =  M33] 

* 

MDIAG  =  0.5*(RPRP**2.0  +  RPAR**2.0) 

MOFFD  =  0.5*(RPRP**2.0  -  RPAR**2.0) 

MDIAQ  =  -RPRP*RPAR 

* 

*  Multiply  reflection  Mueller  matrix  by  Stokes  vector 

* 

ITMP  =  MDIAG*ISRF(I,J)  +  MOFFD*QSRF(l,J) 
QTMP  =  MOFFD*ISRF(I,J)  +  MDIAG*QSRF(I.J) 
UTMP  =  MDIAQ*USRF(I,J) 

* 


*  Specify  Mueller  matrix  for  depolarizing  interface 

*  [MDIAG  =  M22  =  M33,  Ml  1  =  1 .0] 

* 

MDIAG  =  DEPOL 

* 

*  Multiply  depolarization  Mueller  matrix  by  Stokes  vector 

* 

ISRF(I,J)  =  ITMP 
QSRF(I,J)  =  MDIAG*QTMP 
USRF(I,J)  =  MDIAG*UTMP 

•k 

*  Rotate  Stokes  vector  back  to  reference  plane  of  polarizer 

* 

QTMP  =  COS(ALPHrQSRF(l,J)+SIN(ALPH)*USRF(l,J) 
UTMP  =  -SIN(ALPHrQSRF(l,J)+COS(ALPHrUSRF{l,J) 
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QSRF(I,J)  =  QTMP 
USRF(I,J)  =  UTMP 

* 

*  Calculate  reflected  surface  intensity  from  reflected  radiance  via  finite  surface 

*  facet  attenuation  and  projection  to  sensor: 

*  ATT(omega,beta)  =  cos(omega)*sec(beta)  =  cos(omega)/cos(beta) 

* 

IF  (BETA.LT.HALFPl)  THEN 

ATT  =  COSOMEGA/COS(BETA) 

ELSE 

WRITE(*  *)  'Beta  >  90  deg:  surface  Is  non-analytic' 
WRITE(*,*)  'setting  radiance  attenuation  coefficient  to  O' 
ATT  =  0.0 

ENDIF 

* 

ISRF(I,J)  =  ISRF(I.J)*ATT 
QSRF(1,J)  =  QSRF(I,J)*ATT 
USRF(I.J)  =  USRF(I,J)*ATT 

* 

PSRF(I,J)  =  SQRT(QSRF(I,J)**2.0  +  USRF(I,J)**2.0)/ISRF(I,J) 
TSRF(I,J)  =  0.5*ATAN2(USRF(I,J).QSRF(I,J)) 

* 

*  Multiply  non-ideal  linear  polarizer  Mueller  matrices  by  Stokes  vector 

*  for  four  polarizing  angles:  0, 45,  90,  and  135  degrees  -  to  obtain 
*l(0),  l(45),  l(90),  and  l(135) 

I000(I,J)  =  0.5*(K_1  +  K_2)*(ISRF(I,J)  +  QSRF(I,J)) 

I045(I,J)  =  0.5*(K_1  +  K_2)*(ISRF(I,J)  +  USRF(I,J)) 

I090(I,J)  =  0.5*(K_1  +  K_2)‘(ISRF(I,J)  -  QSRF(I,J)) 

I135(I,J)  =  0.5*(K_1  +  K_2)*(ISRF(I,J)  -  USRF(I,J)) 

* 

*  Calculate  maximum  intensities  for  each  of  the  intensity  distributions 
*1(0),  1(45),  1(90),  and  1(135) 

*  [l(0)  &  l(90)  are  Image  pairs;  l(45)  &  I9135)  are  image  pairs] 

* 

IF  (I000(I,J).GT.IMAX1)  IMAX1  =  I000(I,J) 

IF  (I090(I,J).GT.IMAX1)  IMAX1  =  I090{I,J) 

IF  (I045(I,J).GT.IMAX2)  IMAX2  =  I045(I,J) 

IF  (I135(I,J).GT.IMAX2)  IMAX2  =  I135(I,J) 

* 

4  ENDDO 

ENDDO 

* 

WRITE(*,*)  • ' 

WRITE(*,*)  ’Starting  2nd  loop  ...  ’ 

WRITEf,*)  ■|MAX1,IMAX2  =  ’,IMAX1,IMAX2 
WRITE(*,*) ' ' 

* 

DO  I  =  P+1, Q 
DO  J  =  P+1,Q 

★ 

*  Calculate  effective  Intensity  measured  by  sensor:  effect  of  film  exposure 

*  +  effect  of  film  processing  +  effect  of  digital  transformation 

*  [Leff  =  DN_max*((l+bias)^(-gamma_eff))] 

* 

I000(I,J)  =  I000(l,j)  +  BIAS 
I045(I,J)  =  I045(I,J)  +  BIAS 
I090(I,J)  =  I090(I,J)  +  BIAS 
I135(I,J)  =  I135(I,J)  +  BIAS 

* 

I000(I,J)  =  FLOAT(INT(DN_MAX’‘((IOOO(I,J)/IMAX1)**(-GAMMA)))) 
I090(I,J)  =  FLOAT(INT(DN_MAX*((I090(I,J)/IMAX1)**(-GAMMA)))) 
1045(1, J)  =  FLOAT(INT(DN_MAX*((I045(I,J)/IMAX2)**(-GAMMA)))) 
I135(I,J)  =  FL0AT(INT(DN_MAX*((I135(I,J)/IMAX2)**(-GAMMA)))) 


*  Recalculate  Stokes  images  from  intensity  images 
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ISRF(I.J)  =  0.5*(I000(I,J)  +  I045(1.J)  +  1090(1,  J)  +  I135(I.J)) 

QSRF(I,J)  =  I000(I,J)  -  I090(I,J) 

USRF(I.J)  =  I045(I,J)  - 1135(1, J) 

* 

IF  (ISRF(l,J).NE.0.0)  THEN 

PSRF(I,J)  =  SQRT(QSRF(I.Jr2.0  +  USRF(I.Jr2.0)/ISRF(l,J) 

ELSE 

PSRF(I,J)  =  0.0 

ENDIF 

* 

*  Recalculate  Stokes  derivative  images  from  Stokes  images 

* 

IF  (USRF(1,J).NE.0.0)  THEN 

TSRF(I,J)  =  0.5*ATAN2(USRF(I,J),QSRF(I,J)) 

ELSE 

TSRF(l,J)  =  0.0 

ENDIF 

* 

ENDDO 

ENDDO 

* 

RETURN 

END 
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